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With rapid progress made in employing computational techniques for various 
complex Navier-Stokes fluid flow problems, design optimization problems traditionally 
based on empirical formulations and experiments are now being addressed with the aid of 
computational fluid dynamics (CFD). To be able to carry out an effective CFD-based 
optimization study, it is essential that the uncertainty and appropriate confidence limits of 
the CFD solutions be quantified over the chosen design space. The present dissertation 
investigates the issues related to code verification, surrogate model-based optimization 
and sensitivity evaluation. 

For Navier-Stokes (NS) CFD code verification a least square extrapolation (LSE) 
method is assessed. This method projects numerically computed NS solutions from 
multiple, coarser base grids onto a finer grid and improves solution accuracy by 
minimizing the residual of the discretized NS equations over the projected grid. In this 
dissertation, the finite volume (FV) formulation is focused on. The interplay between the 
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concepts and the outcome of LSE, and the effects of solution gradients and singularities, 
nonlinear physics, and coupling of flow variables on the effectiveness of LSE are 
investigated. 

A CFD-based design optimization of a single element liquid rocket injector is 
conducted with surrogate models developed using response surface methodology (RSM) 
based on CFD solutions. The computational model consists of the NS equations, finite 
rate chemistry, and the k- s turbulence closure. With the aid of these surrogate models, 
sensitivity and trade-off analyses are carried out for the injector design whose geometry 
(hydrogen flow angle, hydrogen and oxygen flow areas and oxygen post tip thickness) is 
optimized to attain desirable goals in performance (combustion length) and 
life/survivability (the maximum temperatures on the oxidizer post tip and injector face 
and a combustion chamber wall temperature). A preliminary multi-objective optimization 
study is carried out using a geometric mean approach. Following this, sensitivity analyses 
with the aid of variance-based non-parametric approach and partial correlation 
coefficients are conducted using data available from surrogate models of the objectives 
and the multi-objective optima to identify the contribution of the design variables to the 
objective variability and to analyze the variability of the design variables and the 
objectives. 

In summary the present dissertation offers insight into an improved coarse to fine 
grid extrapolation technique for Navier-Stokes computations and also suggests tools for a 
designer to conduct design optimization study and related sensitivity analyses for a given 
design problem. 
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CHAPTER 1 

INTRODUCTION AND SCOPE 

Introduction 

With the advancements in computational technologies, different complex design 
problems can now be analyzed with reduced economical and computational cost as 
compared to an experiment. Additionally, computational models can be used as a 
substitute to study environments that are too hazardous to conduct physical experiments. 
These aspects have motivated the development of different tools in the areas of 
computational fluid dynamics (CFD) and multi-disciplinary optimization (MDO) which 
in turn promotes CFD-based design optimization studies. 

Computational fluid dynamics-based design studies can be broadly divided into 
two parts. The first part is the accurate modeling of the fluid flow problem and the second 
is the sensitivity estimation of the CFD solutions to the design variables and exploration 
of optimum designs. Solving fluid flow problems based on Navier-Stokes (NS) equations 
and other related models introduces challenges that include physical modeling 
uncertainty, geometric complexities, non-uniform and non-orthogonal meshing, disparate 
length scales and characteristics of the flow variables (such as velocity and pressure), and 
acceptable run time for engineering analysis and design. To obtain an accurate numerical 
solution after addressing all these issues is an elaborate topic of research. As pointed out 
in the review by Oberkampf and Trucano (2002), there are 5 sources of error in a CFD 
analysis, namely, insufficient spatial discretization convergence, insufficient temporal 
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discretization convergence, insufficient convergence of an iterative procedure, computer 
round-off and computer programming error. Much work has been done to estimate the 
first three errors. The remaining two are related to the field of computer programming 
and usually involve the verification of the code or the software. 

In order to ascertain the overall accuracy and identify uncertainties associated with 
physical modeling there is a strong need to develop techniques for accuracy assessment 
of a given numerical approach and mesh resolution. A priori and a posteriori error 
estimates are widely used to estimate errors and verify CFD solutions. A priori error 
estimation uses the information available before a solution is obtained based on the 
partial differential equation (PDE) that has to be solved and the initial and boundary 
conditions. In CFD analysis these estimates aid in the estimation of the stability and 
existence of a solution. On the other hand, an a posteriori error estimate uses the 
information of the computed solution. This allows the performance of a given numerical 
approach to be judged by estimating the errors arising out of the discretization of the 
PDEs. 

The issue of mesh resolution is of importance because a coarse grid is 
computationally economical but maybe inadequate in capturing the detail flow features 
whereas a fine grid, although it provides a well-resolved solution, increases the 
computational cost. Hence it is of interest to use a coarse grid and obtain an accurate 
solution. A well known method in this category is the Richardson extrapolation (RE). 
Details of this method can be found in Roache (1997) and Shyy et al. (2002). The 
Richardson extrapolation is based on eliminating the leading order error term based on 
the local truncation error analysis, implemented via solutions obtained on two or more 
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base grid levels. An attractive feature of this method is that it can be applied to any partial 
differential equation (PDE) without direct information about the solution procedure or the 
equation itself. On the other hand, the disadvantage of this method is that it assumes the 
order of the solution a priori, whereas generally the order of convergence of a CFD 
solution is space and solution profile dependent, and may not be consistent with that 
indicated by the local error analysis. As pointed out by Shyy et al. (2002), the flow field 
with high gradient in flow properties can deteriorate the performance of this method. 
Typically, RE is usually conducted on two base grids with one having twice the number 
of nodes as the other in each direction. For practical problems with complex geometries, 
it is often difficult to satisfy such a requirement. A very coarse grid may not efficiently 
capture the flow features for extrapolation and a very fine grid might increase the 
computational cost, thereby spoiling the whole purpose of extrapolation. 

Recently, Garbey and Shyy (2003) have proposed a least square extrapolation 
(LSE) technique. This method addresses some of the disadvantages noticed in RE. Unlike 
RE, LSE directly uses the discretized PDE model under consideration to estimate the 
degree of convergence. In this method residual is estimated on the fine grid onto which 
the base grid solutions are extrapolated based on the discretized equation and the specific 
computational scheme. Additionally the space dependency of the solution’s convergence 
order can be accounted for by assuming spatially dependent weight functions. Garbey 
and Shyy (2003) have shown that LSE is effective even when the base grids are not based 
on grid doubling or halving. Details of the LSE method will be discussed in chapter 2. 

In the area of MDO tools like design of experiment (DOE) (Myers and 
Montgomery, 1995) and response surface methodology (RSM) (Myers and Montgomery, 
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1995) have been shown to work well for a wide range of design problems. Hence it is of 
interest to use these tools to design devices that have complex geometries and flow 
features. Sensitivity analysis and optimization studies related to these devices require 
numerous function evaluations and the computational cost of CFD simulations restricts 
its use in such scenario. For such situation, RSM has been found to be very useful. A set 
of designs can be identified using DOE techniques and the solution for the designs 
obtained using CFD analyses. These solutions are then represented with response surface 
approximations (RS As) (polynomials or any other analytical function), which can then be 
used as the function evaluator. 

Sensitivity analysis helps the designer by providing the knowledge of the influence 
of design variables on the essential flow features. Global sensitivity analysis (Sobol, 
1993) is a method which has the potential for providing such information. The global 
sensitivity indices can be used to rank the design variables in order of their importance 
and also to identify the interaction between them. Based on these observations, 
nonessential variables can be fixed and the dimensionality of the design problem reduced. 

Most of the design studies require that more than one aspect of the solution 
(objective) be improved upon at the same time. There are different methods available for 
solving such multi-objective problems. Multi-objective optimization problems usually 
have numerous optimal solutions, known as Pareto optimal solutions (Miettinen, 1999). 
These solutions are such that no improvement is possible in any objective without 
sacrificing at least one of the remaining objectives. Hence these solutions are non- 
dominated. On the other hand if for a given solution there are other solutions where 
improvement in any objective is possible without sacrificing the remaining objectives, 
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that solution is said to be dominated. The function space of all the non-dominated 
solutions in the Pareto-optimal set is termed the Pareto-optimal front (POF). This set 
provides the designer with a clear idea of the trade-offs involved in a design study. 
Evolutionary algorithms (EAs) are popular global optimizers that have been used to find 
multiple Pareto optimal solutions in a single simulation run. Different multi-objective 
evolutionary algorithms (MOEA) have been proposed in literature (Deb, 2001). Insight 
into the trade-offs between different objectives can be obtained from the POF. 

In this dissertation different aspects of NS code verification and CFD-based design 
process are addressed and methods explored with the help of individual case studies. The 
LSE used for NS code verification is implemented on a lid-driven cavity flow with 
different boundary conditions and Reynolds numbers. The goal is to assess the 
effectiveness of the LSE technique for finite volume (FV) NS computations. 

The CFD-MDO coupled design study is addressed with the aid of a single element 
liquid rocket injector. The study is conducted using response surface methodology along 
with global sensitivity analyses and genetic algorithms. The overall goal is to estimate the 
global sensitivity and trade-offs involved in the design of a rocket injector. 

Scope 

The scope of the dissertation can be broadly divided into two, namely, 

• Investigation of NS CFD code verification. 

• CFD-based design optimization. 

The NS CFD code verification deals with the demonstration of LSE for numerical 
accuracy improvement and computational cost reduction with the aid of few widely used 
benchmark problems. The implementation of the method for complex flows involved in 
injector design is an issue of future study. The goal of the LSE related study is to 



6 


• Assess the effectiveness of the LSE technique for FV computations. 

• Evaluate the implication of solution gradients and singularities on the performance 
of LSE. 

• Address the issue of coupling between the pressure and velocity in the NS 
equations. 

The CFD-based design optimization study is addressed by modeling a single 
element liquid rocket injector. The goal is to integrate CFD with the optimization tools 
and obtain better design methodologies than those used in the past. While the ultimate 
goal is to analyze multi-element injectors, much of the detailed work in injector design 
can be done, or at least initiated, at the single element level. 

Design variables governing the life/survivability and performance of the injectors 
are identified, namely, H 2 flow angle, H 2 and 0 2 flow areas with fixed mass flow rates of 
fuel and oxidizer and 0 2 post tip thickness. The design objectives that are 
life/survivability indicators are the maximum temperature on the oxidizer post tip, the 
maximum temperature on the injector face and the combustion chamber wall temperature 
taken three inches from the injector face. The performance indicator is the length of the 
combustion zone. 

To facilitate the development of the present methodology, a baseline element 
design is needed. This baseline concept is generated by an empirical design methodology 
based on a specific set of propellant flow rates, mixture ratio and chamber pressure. The 
selected design variables are then varied based on this baseline design and the design 
space populated with the aid of a DOE technique. The prescribed CFD cases are executed 
and post processed to extract the required dependent variable data. This data are then 
used to generate a RSA for each objective in terms of the independent design variables. 
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Sensitivity and trade-off analyses for the design are then carried out using the data 
obtained from surrogate models of the design objectives, and the POF (multi-objective 
optima) generated by Goel et al. (2004) with the aid of a multi-objective genetic 
algorithm (MOGA) and a local search method (e -constraint strategy). The regions of the 
POF that represent different trade-offs among the objectives are obtained through a 
hierarchical clustering algorithm which will be discussed in chapter 3. This study is 
broadly divided into the following two parts. 

1. A sensitivity study is carried over the whole design space. The contribution 
of the design variables to the objectives variability is calculated using a 
variance-based non-parametric approach and correlations between objectives 
are investigated. 

2. Sensitivity analyses are conducted on clusters along the POF. Box plots are 
used to highlight the variability of the design variables and the objectives 
within a cluster. Additionally, the linear relationships between the design 
variables and the objectives are explored with the aid of partial correlation 
coefficients. 

The remaining chapters of this dissertation address each of the issues mentioned 
above in detail. In chapter 2 the different aspects of NS CFD code verification are given. 
In this chapter a literature survey is presented initially followed by the details of 
Richardson extrapolation (RE) and least square extrapolation (LSE) techniques. These 
techniques are implemented on a two-dimensional turning point problem and a laminar 
lid-driven cavity flow. The performances of the RE and LSE are evaluated and issues 
related to their implementation are identified. 

In chapter 3 different tools used for the design optimization study are presented in 
detail. Response surface methodology, DOE techniques, MOGA, optimization algorithm 
and sensitivity analyses tools are presented. Additionally a clustering method and a 
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visualization tool are presented which along with other tools mentioned above are used in 
chapter 4 to study the design problem in detail. 

In chapter 4 the past and present studies related to the design of a single element 
rocket injector are presented. The design variables and objectives of the current injector 
design are identified and optimization studies carried out for different design goals using 
surrogate models developed based on CFD solutions. Sensitivity analyses are done with 
the aid of surrogate models and sensitivity indices. The design trends are observed with 
the aid of results obtained from these studies and visualization techniques. 

In chapter 5 the dissertation is summarized and conclusions drawn. The unresolved 
issues are identified and scopes for future studies identified. 



CHAPTER 2 

NAVIER-STOKES CODE VERIFICATION 
In this chapter, the literature survey related to the work done on computational fluid 
dynamics (CFD) code verification is presented. Following this, details of Richardson 
extrapolation (RE) and least square extrapolation (LSE) techniques are presented along 
with case studies and results. 

Literature Review 

In the area of CFD code verification, there has been considerable research done to 
address the issues of grid sensitivity and quantification of uncertainties in the 
computations (Roache, 1997; Pelletier et al., 2001; Gu et al., 2001; Turgeon et al., 2001). 
The governing partial differential equations (PDEs) of fluid mechanics are either linear or 
nonlinear depending on the problem solved. The most frequently used computational 
methods for solving these PDEs are finite element (FE), finite difference (FD) and finite 
volume (FV) methods. To obtain an accurate solution using any of these discretization 
schemes, adequate grid resolution is essential since for a consistent scheme the computed 
solution approaches the exact continuous solution as the grid spacing tends to zero (Shyy, 
1994). The error arising out of inadequate grid distribution falls under the category of 
domain discretization error. At the same time the discretization operators used for the 
derivatives in the PDEs should be according to the order of variation of the solution over 
the domain. If an improper variation is assumed, then it gives rise to the equation 
discretization error. Hence the goal is to achieve a fine grid distribution over the domain 
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and a valid discretization of the derivatives in the PDEs without violating issues related to 
the stability of the computational model. Additionally, to obtain a solution for these PDEs 
valid boundary and initial conditions need to be specified which will define the flow 
field. 

The error estimation can be done before or after the CFD analysis is carried out. 

The former is called a priori error estimate and the later is referred to as a posteriori error 
estimate. In this dissertation the aim is to implement the a posteriori error estimates to a 
Navier-Stokes (NS) CFD code. Hence a literature survey of the available a posteriori 
error estimates is presented. As pointed out by Oden (2001), there are two broad classes 
of error estimation methods, namely, residual methods and recovery methods. In the 
residual method, the residual estimates the degree with which the approximate solution 
fails to satisfy the equations of the original problem (Babuska and Rheinbolt, 1978; 
Demkowicz et al., 1984). The residual estimates are used to define error norms which are 
largely used for controlling adaptive meshing procedures. The other type, referred to as 
the recovery method, attempts to enhance the computed solution during the post- 
processing step, like Richardson extrapolation (Roache, 1997; Shyy et ah, 2002) and least 
square extrapolation (Garbey and Shyy, 2003). 

There is a vast amount of literature available with respect to a posteriori error 
estimators for finite element approximation. A recent book by Ainsworth and Oden 
(2000) provides a survey of the available literature on such error estimators. Babuska and 
Rheinbolt (1978), Demkowicz et ah (1984), Zienkiewicz and Zhu (1987) are some of the 
earliest works in this area. These works concentrate on error estimation in finite element 
problems. Zienkiewicz and Zhu (1992a, 1992b) in a two-part paper introduced a recovery 



11 


technique which was an improvement on their previous effort. They achieved higher 
order accuracy by using this approach and also found it to be cost effective. Their 
approach uses a single and continuous polynomial expansion for the nodal values of the 
solutions (displacements, stresses) over a patch of elements surrounding the node of 
interest. A least square method or an L 2 projection is used to obtain this polynomial. The 
error in the computed solution is then measured by comparing it with the recovered 
solution using different forms of error norms. Ainsworth and Oden (1992, 1993a, 1993b) 
have presented an element residual method for finite element computations. This method 
has later been implemented by Jasak and Gosman (2001) for finite volume approach. 
Oden et al. (1994) extended the application of element residual methods to NS equations 
using a finite element approach. They designed error norms to measure the error in 
velocity and pressure for incompressible flow problems. 

Use of a posteriori estimator in finite volume approach has been of interest only 
over the past few years. Ilinca et al. (2000) have compared three different error estimators 
for finite volume solutions. They have compared a technique using Richardson 
extrapolation, a technique based on the difference of the computed solution, and a higher- 
order reconstruction obtained using a least square method and a technique which solves 
an error equation. Richardson extrapolation (Roache, 1994) uses the solutions on 
different grids with different levels of refinement, one finer than the other. Using the 
Taylor series representation of the discrete solution and combining the available multiple 
solutions, a higher order approximation of the desired variable is obtained. The drawback 
is to generate an adequately resolved solution on multiple grids for complex problems 
and the a priori assumption of the solution accuracy order. Solution reconstruction 
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(Barth, 1993; Ollivier-gooch, 1996) based on an approximation of the derivatives in the 
Taylor series expansion using a weighted least square method improves solution 
accuracy. This is then used to compare with the initial solution to obtain an error 
measure. The third method proposed by Zhang et al. (1997) estimates error using an error 
equation. In this method the error equation is derived directly from the governing 
equation and has a source term, which needs to be computed. This source term is 
evaluated using the flux differences at the cell interface. These methods were tested on 
subsonic, transonic and supersonic flows. Richardson extrapolation and the error equation 
method are shown to perform reasonably well as compared to the solution reconstruction 
method. 

Richardson extrapolation is based on eliminating the leading order error term in the 
assumed Taylor series expansion of the solution. This is done by combining solutions 
obtained on two different grids with (uniform) discrete spacings. The attractive feature 
about this method is that it can be implemented on a solution for any PDE without 
bothering about the solution procedure or the equation itself. On the other hand, the 
disadvantage of this method is based on the fact that it assumes a fixed order of accuracy 
(when using two grids) of the solution all over the computational domain. In practical 
fluid problems, due to numerical schemes and fluid physics involved, the order of 
accuracy is not uniform over the computational domain. Hence, the local truncation error 
based on Taylor series expansion may not represent the global accuracy of the solution. 
To address this, three grids of different resolution are required to estimate the order of 
accuracy of the solution. Additionally, problems involving turbulent flow features will 
result in sharp gradients in flow parameters in different regions of the domain. It has been 
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noticed by Shyy et al. (2002) that flow field with high gradient in flow properties can 
deteriorate the performance of this method. Although the most common use of RE is with 
two grids, one twice as fine as the other, this is not always essential. But for coarse grids, 
the assumption of monotonic truncation error convergence may not be valid thereby 
requiring fine base grids. For practical problems with complex geometries, it is tough to 
obtain two grids that satisfy such requirements. A very coarse grid may not efficiently 
capture the flow features for extrapolation, and a very fine grid might increase the 
computational cost, thereby spoiling the whole purpose of extrapolating. Another 
important disadvantage of RE is that the extrapolation does not preserve the 
conservativeness of the flow properties. 

The book by Roache (1998) has more information on different approaches to 
address some of these issues. But in its simplest form, RE is based on an a priori 
assumed order of convergence of the continuous solution. In the work of Garbey and 
Shyy (2003), basic properties of RE and its computational implications are presented in 
detail. They have presented complementary views that are asymptotic expansion for 
continuous function in a normed vector space and numerical approximation for discrete 
functions defined in a mesh. Comparing these views they have pointed out that it is 
difficult to achieve asymptotic order of convergence unless the numerical perturbations 
(arising out of imperfect convergence of the fluid solver) in the computations are small. 
Unless a grid is very fine, it is hard to reduce the perturbation error. Based on their study 
on the convergence order approximation and RE in relation to CFD problems, they have 
found that convergence order is space and grid dependent. Hence, they have concluded 
that there is no way to justify the efficiency of RE which uses a constant weight 
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formulation. Additionally, RE is easy to implement only when the grid spacing is 
uniform and the subsequent grids used in the extrapolation are obtained by grid doubling 
or halving. 

Scope 

The aim of this work is to investigate the performance of LSE method in NS 
computations (Vaidyanathan et al., 2004a). In particular, FV formulation is considered, 
which is different from the FD approach adopted by Garbey and Shyy (2003). The lid- 
driven cavity flow with different boundary conditions and Reynolds numbers (Re) is 
adopted as the main test problem. The goal of this study is to address the following 
issues: 

• To assess the effectiveness of the LSE technique for FV computations. 

• To evaluate the implication of solution gradients and singularities on the 
performance of LSE. 

• To address the issue of coupling between the pressure and velocity in the NS 
equations. 

It is noted that the geometric complexity plays a major role in practical flow 
computations. This aspect can have substantial impact on the performance of any 
extrapolation techniques. This issue will be addressed in future studies. 

In the following sections, the salient features of both RE and LSE techniques will 
be first presented. Following this, a brief discussion of the flow solver and the algorithms 
for the implementation of the LSE in relation to the NS computations will be discussed. 
Finally, their implementation on test cases and a discussion based on the results obtained 
will be provided. 
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Approach 

In the following, RE is reviewed first to help motivate the main topic of interest, 
namely the least square extrapolation. 

Richardson Extrapolation (RE) 

In essence, RE is based on an a priori assumed order of accuracy of the continuous 
solution. Considering a flow-field output, say a velocity component or pressure, of a CFD 
simulation on a grid size, h, and assuming second-order convergence a priori , the Taylor 
series expansion can be written as 

U(x) = u(x;h) + a 2 h 2 +a 3 h 2 +... (2.1) 


where U(x) is the exact solution and u(x;h) is the numerical solution based on h. Similarly 
the solution on a grid size, h/2 can be written as 


U{x)-u 


h\ 

x ’2 

l ) 


h 2 h 

+ a 2-J + C1 lJ + - 


( 2 . 2 ) 


Solving for U(x), using Eqs. (2.1) and (2.2) and eliminating 0(h 2 ) term and 
neglecting the higher-order terms leads to 
1 


U(x)=- 


4u 


+ O (A 3 ) 


(2.3) 


This is a second-order RE which results in a gain of order one. Similarly if the 
solution is assumed to have first-order convergence a priori and 0(h) term is eliminated, 
it results in a first-order RE with a gain of order one. 


U(x) = 2u x;^-u(x;h) + o(h 2 } 


(2.4) 


When in a practical fluids problem, the degree of convergence is not uniform over 
the computational domain, the order of extrapolation used can either improve the order of 
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accuracy or not affect it at all. If Eq. (2.4) was based on assumptions made in Eqs. (2.1) 
and (2.2), there would be no gain in the degree of convergence. Equations (2.3) and (2.4) 
can be generalized for two grids with spacing hi and h 2 , resulting in the following 
respective extrapolation scheme. 

u(x)= !^MtM!L S. ( 2 . 5 ) 


hfi(x,h,}-h,u(x,h.) _ 

u(x) = * K — y - ( 2 - 6 ) 

\K~h 2 ) 

where u(x,hj) represents the interpolated values from two coarse grid solutions, which are 
not necessarily based on grid doubling or halving. In implementing the extrapolation 
techniques, it is required that the coarse grid solutions be interpolated to the locations of 
the fine grid on which the extrapolated solution is obtained. The order of interpolation 
has to be such that it does not deteriorate the order of the computed solution. 

Least Square Extrapolation (LSE) 

In this dissertation the least square extrapolation method is used which, most 
importantly, estimates the order of convergence as the solution to a least square problem 
instead of assuming it a priori. Additionally, it accounts for the dependence of the order 
of accuracy on spatial coordinates by using spatially dependent extrapolating weights. 
The details of this method are presented below. In this approach, two coarse grid 
solutions, ui and U 2 , not necessarily based on grid doubling or halving, are combined 
linearly using a weighting function, a, which can be spatially dependent. The 
extrapolated solution is given as 

U =au l +(l-a)u 2 (2.7) 
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where u/ and u 2 represent the coarse grid solutions interpolated onto a fine grid, M 0 . 

Following Garbey and Shyy (2003), a modified Fourier expansion has been used 
for the weighting function, a. The weighting function, a is given by 

M 

a (jc) = a 0 + a x cos(x;r) + sin {(j- 1 (2.8) 

with (Xj,j = reals. The weighting function, oris dependent on the spatial co- 

ordinate, x for a one-dimensional problem. For two dimensional problem it is a function 
of spatial co-ordinates, x and y such that 

a(x,y) = a(x)a(y) (2.9) 

This modified Fourier expansion is ideal for rectangular domains. Different function may 
have to be used for more complex domains. The coefficients, ay are estimated by solving 
a least squares problem. For a given linear PDE, say 

L[u] = f (2.10) 

with u e (£ a ,|| || a ) and / € [E b ,\\ || fc j , its numerical approximation can be written as 

H=/* 

with u h e (£*,|| and f h e (£*,|| ||J . Based on this, the least square problem can be 
formulated as 

P a : Find a e /1(0,1) c L x such that 

L h [au x +(l-a)u 2 ]-f h (2.12) 

is minimum in L 2 (M C ,), where M 0 represents the fine grid to which the coarse grid 
solutions are interpolated. For a one-dimensional problem using two Fourier modes, Eq. 
(2.12) that has to be minimized can be rewritten as 
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L h [a 0 «, + (1 - a 0 )u 2 ] + . L h [a, cos {xn)u x + (l - a, cos(x;r ))u 2 ] ■ - /* (2.13) 

This approach can be generalized for nonlinear PDE problems via a Newton-like loop. 
For a given nonlinear problem, say, 

W[«] = / (2-14) 

linearization results in 

/(«)[v] = g (2-15) 

The least square problem then becomes 

J h {a k u l + (l-a*)w 2 )[a* +l M, +(l -a* +1 )« 2 ]-g A (2.16) 

which is minimum in L 2 (M 0 ). The iteration is started from the initial condition of = 0, 
until ||q:* +1 -a k | is less than some tolerance value. This method requires that the initial 

guess should not be too far from the final solution for convergence; i.e., the solution ui 
should be close to the grid solution on M 0 . In situations where J(u) is not available, a 
nonlinear set in a is obtained, which can be solved using nonlinear solvers. This is an 
issue of further study. In this initial implementation, a linearization is adapted to simplify 
the problem in hand. 

A point to be noted is that LSE uses interpolated functions on a fine grid, and the 
differential order, say n, of the PDE that is being solved impacts the choice of the 
interpolating scheme. The scheme used should give a result that should be smooth 
enough to be n-times differentiable. In the present study, cubic spline is used for 
interpolation of the solutions. 

The LSE technique is based on the minimization of the residual. However, the real 
goal of LSE is to maximize the solution accuracy, or equivalently, minimize the solution 
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error on the extrapolated grid. However, minimization of the residual is not necessarily 
the same as the minimization of the solution error. This aspect is explained in the 
following subsection. 

Error versus Residual 

To see the above mentioned point, let us consider only the linear problem, 
assuming the following form: 

AU = b (2.17) 

where A is a symmetric coefficient matrix of size nxn and U is the n x 1 exact solution 
vector of the linear system. Assuming u to be a non-exact solution of size n x 1, we have 
the following relationship: 

Au = b + r (2.18) 

where r represents the n x 1 residual vector of the system of equations. Defining e as the 
difference between u and U and subtracting Eq. (2.18) from Eq. (2.17) gives: 

As = -r (2.19) 

Denoting the eigenvectors of matrix A by a„ i = 1, ..., n such that ||a,|| = 1.0 and the 
eigenvalues as i = 1, ..., n the following relationships can be obtained: 

r = ( 2 - 20 ) 

/=1 

where i= l, ...,n are linear weights associated with each eigenvector and 

e = -A'‘r (2.21) 

which can be rewritten as 

<s = E(«,Mk (2.22) 

i-i 


The residual, F, is then defined as 



The least square approach aims at minimizing the I 2 norm of the residual, F; i.e., 


( n 

min (F) = min \ ^ |/; | 


( n 


= mm 


Ehf 


(2.24) 


v 1=1 


Now the adequacy of the non-exact solution, u, can be measure by the L 2 norm 
error measured between the u and U, which can be written as: 


error = E 


n 

=2>, 


(2.25) 


1=1 


Minimizing the error, gives 


n 

min (error) = min ^ |(a,. / )| 

l /*=! 


2A 


(2.26) 


From Eqs. (2.24) and (2.26) it is clear that minimization of F is equivalent to 
minimization of L 2 error norm {error) only when A,-’s are equal. Therefore it is clear that 
properties of the matrix A define the properties of 2.,’s and the relationship between F and 
L 2 error norm {error). The implication of the residual minimization in accuracy control 
will be examined along with the test problems. 

In the context of the current study, the solution error measure is defined as follows 


e: = 




1/2 


j= 1 1 

where n indicates the source of the non-exact solution, m indicates the source of the 
solution with respect to which the error is measured; i.e., the reference grid, (u.j ) 


(2.27) 


represents the extrapolated or interpolated solution (non-exact solution) onto the fine 
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grid, M 0 , at the (i,j) th node, ) represents the reference grid solution, at the (i,j)‘ h node 


and N is the number of nodes of the fine grid (M 0 ). The reference grid can be either the 
fine grid ( M 0 ) solution or a fixed grid solution (M„) fine enough to be considered as the 
benchmark solution. 

Flow Solver 

The governing equations are solved using the in-house code STREAM which is 
based on the pressure-based algorithm and finite volume approach (Thakur et al., 2002). 
The equations are solved on a multiblock structured curvilinear grid. The convective 
terms in the momentum equations are discretized using the second-order upwind scheme 
(Shyy, 1994) and the diffusion terms are discretized using a central difference scheme. 
The continuity and momentum equations in Cartesian coordinates are presented below. 
Steady state computations are carried out. 

Continuity: 

T-(/»)+T-(/»)-° < 2 - 28 > 

OX ox 

Momentum: 


— (pu(j)) + — (pv<f>) = — 
dx y ' dy KH ’ dx 


f 


df 


d<f> 


^ dx) + dyy dy j 


-S D 


(2.29) 


where p is the density, u and v are the velocity components in Cartesian coordinates, p is 
the laminar viscosity and Sp represents the pressure gradient term. In Eq. (2.29), <j> 
represents u or v component of the velocity. 


Figure 2.1 shows the collocated grid system. All the flow characteristics, namely 
pressure, velocity, density and viscosity are stored at the cell center (P). Carrying out the 
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volume integral on Eq. (2.29) and using the notations shown in Figure 2.1 for the nodal 
references, the following discretized equation is obtained for the momentum equations. 
Gptfip — &E$E a J.W — Sp + S c (2.30) 

where the a ’s represent the convective and diffusive fluxes at the respective grid 
locations, Sp is the source term arising out of the pressure gradient term and Sc is the 
source term arising out of higher order derivatives in the convective fluxes. More details 
of this equation can be found in Shyy (1994) and Shyy et al. (1997). The residual in the 
context of LSE is defined by subtracting the LHS from the RHS of Eq. (2.30). The flow 
solver is based on the SIMPLEC (Van Doormal and Raithby, 1984) algorithm where the 
pressure is updated by solving a pressure correction equation formulated out of 
combining momentum and continuity equations. For the current study, uniformly spaced 
Cartesian mesh is used since the case studies involve a square domain. 
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Figure 2.1 : Collocated grid and notation for 2D grid. 

In the implementation of LSE, a distinct equation is needed for each of the flow 
parameters that have to be extrapolated. Hence for pressure, the pressure equation is used 
which is obtained by substituting the velocity components, obtained from the momentum 
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equations, into the discretized (in the FV sense) continuity equation. The obtained 
pressure equation is given as 

A pPp ~ A E p E + A w p w + A N p N + A s p s + S (2.3 1) 

where the S is a function of the velocity components and the A’s are obtained by 
modifying the momentum and continuity equations (note that A ’s are different than the 
a’s in the momentum equation). The details of the pressure equation can be found in 
Patankar (1980). In summary, there are three equations to estimate the residuals required 
for LSE, namely Eq. (2.30) for u and v-component of the velocity and Eq. (2.3 1) for the 
pressure. The residuals in the context of LSE are defined by subtracting the LHS from the 
RHS of Eqs. (2.30) and (2.31). 

Algorithm for LSE 

As already mentioned there are issues that need to be addressed while 
implementing the LSE technique, namely, the non-linearity in the momentum equations 
(a’s have velocity components in them which make the equation non-linear) and the 
coupling of pressure-velocity so as to preserve the conservativeness of the flow properties 
in each cell. The pressure equation is linear in nature and therefore is simpler to work on 
provided that the source term is adequately represented. To proceed systematically, the 
Navier-Stokes problem is tackled in two steps. In the first step, only the pressure is 
extrapolated and in the second step the pressure-velocity coupled extrapolation is carried 
out. The algorithms are presented in the following sections. 

LSE for pressure equation only 

This step is to test the efficiency of the approach in extrapolating a given solution 
accurately; therefore the velocities obtained from the CFD solution on the finer of the two 
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coarse grids used for extrapolation are used for the source term in Eq. (2.3 1). The 
extrapolated pressure, plse, is defined as 

Plse ~ a Pi "*■ (l — ® ) Pi (2.32) 

where p { and p 2 are the interpolated pressure from the coarse grids to the fine grid using 
spline. The residual for the LSE method is then obtained by using the pressure equation, 
(Eq. (2.31)). For the preliminary analysis of the Navier-Stokes problem, only one mode 
(Eq. (2.8)) is used to obtain the pressure field. 

The algorithm for the LSE of pressure can be stated as 

• Define p LSE =a p ] + (1 -a)p 2 . 

• Input plse, U 2 > v 2 into Eq. (2.31). 

• Find a that minimizes the residual of Eq. (2.3 1). 

• Calculate plse- 

LSE for coupled pressure-velocity 

The momentum equation (Eq. (2.30)) is nonlinear in nature as the coefficients 
contain the velocity components as well. Hence, to linearize the system, a Picard-like 
iterative scheme (Ferziger and Peric, 1999) is adapted. The equation used for the velocity 
extrapolations comes from Eq. (2.30). 


a n d> n + l =a n f +l +a"f + ' +a n d> n+l +a n f +l -S P + S 

P ' f*LSE E E LSE W *LSE N T N LSE S T Slse ^ 


»+l 

C 


(2.33) 


where Sc +1 is defined based on the new velocity components. 

Now the velocity components and pressure are extrapolated sequentially. Again, 
only one mode (Eq. (2.8)) is used for the extrapolation of all the three parameters. The 
algorithm of the approach is given below 


START 
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• Define plse =ap l + (1 -a) p 2 . 

• Input plse, & 2 , v 2 into Eq. (2.3 1) {ulse and vlse are used from the 2 nd loop). 

• Find a that minimizes the residual of Eq. (2.31). 

• Calculate plse- 

• Define u L se =/ui + (1~Y) &2- 

• Input plse, ulse, v 2 into Eq. (2.33) with (f> replaced by u ( v L se is used from the 2 nd 
loop). 

• Find ythat minimizes the residual. 

• Calculate ulse- 

• Define vlse =P v, + (1-P)v 2 . 

• Input plse, ulse, v L se into Eq. (2.33) with <f> replaced by v. 

• Find p that minimizes the residual of the equation. 

• Calculate vlse - 

• Go back to START for the next iteration. Continue the loop imtil a, y and P 
converge within the specified tolerances. 

These algorithms are so designed that a direct access to the coding of the CFD 
software is necessary. The governing equations have to be modified to accommodate the 
estimation of the weights. 

Results and Discussion 

The least square extrapolation method is tested by implementing it on the following 

cases: 

• Linear, two-dimensional turning point problem (Figure 2.2) and 

• Laminar, lid-driven cavity flow (Figure 2.5). 

1 . variable lid- velocity 

2. constant lid- velocity 
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Two-Dimensional Turning Point Problem 

The linear equation that is solved is given below. 


&14 2 

eAu+a(x,y ) — = 0,xe(0,7r) 
dx 


(2.34) 


where £=0.1 and 


a(x,y) = x- 


n 


+ 0.3 


( 

y - -t 

V 2 y 


(2.35) 


H=0 


u=y(n-y) 

n 


u=-y(7T-y) 


Figure 2.2: 2D turning point problem (dimensions are n x x). 

The boundary conditions of the problem are defined such thatxj; e [0,7i]; u =y(n- 
y) at x = 0 and u = -y(n-y) at x = 7i; u = 0 aty = 0 and n (Figure 2.2). There is no velocity 
component along the y direction. This case has been tested by Garbey and Shyy (2003) 
where they have used a finite difference approach. In this work a finite volume approach 
is used. The problem is defined so that the transition layer of thickness e (where the 
velocity changes direction) centered on the curve a(x,y) = 0 is not parallel to x or y axis 
thereby making the problem two-dimensional. A sample solution is shown in Figure 2.3. 
In the finite volume implementation, central difference is used for the diffusion term and 
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first-order upwind is used for the convection term. Two modest base grids Mi and M 2 of 
sizes 23 x 23 and 29 x 29, respectively are chosen such that there is at least one or two 
grid points in the transition layer for the solution on grid M 2 . Both first-order and second- 
order RE are evaluated and compared with the performance of LSE. For this case study 
LSE is implemented with 4 Fourier modes for a in each spatial variable (instead of one 
mode as in the case of Navier-Stokes problem). The extrapolation is done onto fine grids, 
M 0 , ranging from size 41 x 41 to 101 x 101 in steps of 10. In Figure 2.4, subscripts, RE1 
and RE2 refers to RE of first- and second-order, respectively. The errors are represented 
as E b a which denotes the error in solution obtained using a certain method (RE or LSE) or 
CFD simulation on a given grid (represented as a) as compared to the solution obtained 
through CFD simulation on grid denoted by b (e.g. is the RMS error in u estimated 

using all the points in the computational domain, of M 0 , as compared to the interpolated 
solution from a finer grid (M„) of spacing h/2). 



Figure 2.3: Solution ( u ) contours for the 2D turning point problem, grid = 81 x 81. 
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Figure 2.4: Application: turning point problem with s = 0.1, x-axis is the number of 
points in each direction for the fine grid on which extrapolation is done. 

From Figure 2.4 it is seen that RE2 gives better results than RE1 for grid with poor 
resolution. This is because the transition layer is under-resolved and the dominant error 
comes from the second order diffusion term. As the grid gets finer, performance of 
is better than RE2 as the first order convective error dominates. For some intermedi ate 
grid (51 x 51), there maybe a cancellation between the convection and diffusion errors 
which might result in large improvements for RE as seen in Figure 2.4. Richardson 
extrapolation fails as the grid M 0 gets finer than N = 61. In all the cases, LSE predict the 
fine grid solution with an error less than that of the fine grid, M 0 , which is an 
approximation of the exact solution for N as large as 101. With the modest grid sizes of 
23 x 23 and 29 x 29, there are only one or two grid points in the transition layer which 
has a thickness of s= 0.1. Still, the LSE gives a solution, which is more accurate th an the 
fine grid solution on which the extrapolation is done resulting in a gain of more than one 


order of accuracy. 
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This problem is relatively easy to solve and is used to confirm the potential of LSE 
for FV implementation. Use of 4 modes showed better performance of LSE than using 
one mode. Increasing the number of modes beyond four provided marginal improvement 
in the extrapolated solution. Least square extrapolation is found effective for FV 
computations. Additionally it shows that LSE works well with flows having transition 
layers. For this case the resolution of the coarse grids is adequate. There are still issues 
like nonlinear equations, coupling of flow parameters, sharp gradients and singularities in 
the flow that needs to be addressed. To do so, lid-driven cavity flow is used and the 
details are given in the following section. 

Laminar Lid-Driven Cavity Flow 

Garbey and Shyy (2003) have addressed the lid-driven cavity flow problem in a 
vorticity-streamfunction formulation with a finite difference implementation. In the 
current study the complete NS equations are solved in 2D for the lid-driven cavity flow 
and the pressure along with u and v-components of the velocity field are extrapolated. A 
laminar incompressible flow computation is carried out in a square cavity of dimensions 
1 x 1 for Reynolds numbers 5.33 and 1000. This problem addresses the motion of the 
fluid in a square container induced due to the lid-velocity um in the positive x-direction 
(Figure 2.5). The remaining walls of the container have the no slip condition. 

Two different scenarios are considered. In the first case, the singularity at the end 
points of the sliding lid of the unit square cavity is avoided by choosing a variable speed 
as follows: 

u nd =-16x 2 (l -xf 


(2.36) 
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where uud is the w-component of the velocity at the lid. The v-component at the lid is zero 
and x varies from 0 to 1. In the second case, the lid velocity is kept constant, uuj = 1.0. 
This results in singularities in w-velocity at the lid comers where there is an abrupt 
change of uud from one to zero. 

U=Uud 

v=0 


u=0 

v=0 


. u=0 

* * v=0 

Figure 2.5: Two dimensional lid-driven cavity flow (dimensions 1 x 1). 

The Reynolds number for the flow is defined based on the mean lid- velocity. The 
w-velocity contours are shown in Figures 2.6 and 2.7 for Reynolds numbers 5.33 and 
1000, respectively. Some difference is noticed in the w-velocity contours at the upper 
comer regions of the cavity for both the Reynolds numbers. The negative variable uud 
results in a w-velocity profile which looks like the mirror image of the positive constant 
uud for Reynolds number 1000. The constant distribution of velocity results in the 
presence of singularity at the comers of the lid for the pressure as noticed from Figures 
2.8B and 2.9B. The magnitude of pressure near the comers for the case with variable um 
is about 10 times more than that for the case with constant uud (Figures 2.8 and 2.9). The 
presence of singularity will be an issue while implementing the LSE technique. In the 
later sections, the effect will be highlighted by means of the obtained results. 


u=0 

v=0 

y > 




0.3 



A B 


Figure 2.7: w-velocity contour for grid 171 x 171 and Re = 1000. A) M/,y = -16x 2 (l-x) 2 . B) 

uud= 1 . 0 . 

The positive and negative values of pressure identifies whether the pressure is 
acting on the lid or away from it, respectively. 
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Figure 2.8: Pressure along the lid for grid 121 x 121 and Re = 5.33. A) uud = -16x 2 (l-xf. 
B) uud = 1.0. 




Figure 2.9: Pressure along the lid for grid 171 x 171 and Re = 1000. A) uud = -16x 2 (l-x) 2 . 
B) uud = 1.0. 


The results for the following different scenarios are presented below: 

• LSE for Pressure only (Re = 5.33 and 1000) 

1. uud is variable. 

2. uud is constant. 

• LSE for Pressure- Velocity (Re = 5.33 and 1000) 
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1. uud is variable. 

LSE for pressure only 

The results are presented in two parts. In part (a) the case with variable lid- velocity 
is presented and in part (b) the case with constant lid-velocity is presented. In both the 
parts the results for Re = 5.33 are presented first followed by the results for Re = 1000. 
For Re = 5.33, the two selected base grids which captures the flow features adequately 
are 61 x 61 and 71 x 71. The extrapolated pressure fields are compared to the solution 
obtained using a fixed fine grid (M„) of size 121 x 121. The extrapolated fine grids (M 0 ) 
vary between sizes 81x81 and 121x121. For Re = 1000, the two coarse grids are 1 1 1 x 
1 1 1 and 121 x 121 and the fixed fine grid (M n ) for comparison is 171 x 171. The 
extrapolated grids (M 0 ) are between 131 x 131 and 171 x 171. Only one Fourier mode 
(Eq. (2.8)) is used for the extrapolation of pressure. The results are compared with results 
obtained from Richardson extrapolation. 

Variable lid-velocity (Re = 5.33). In Figure 2.10A the comparison of extrapolated 
solutions on grid M 0 (81 x 81, ...) with the CFD solutions on the same grid (M 0 ) show 
that LSE approach performs better as expected. The error increases as the grid gets finer 
since resolution of the base grid solutions is not adequate for extrapolating onto very fine 
grids. From Figure 2.1 0B, the true error can be estimated where the extrapolated 
solutions on any grid M 0 is compared with the fixed fine grid ( M „ ) solution which 
represents the “correct” solution, obtained on a substantially finer grid, for the given 
Reynolds number. As expected it is seen that as the grid gets finer, LSE outperforms 
other methods and attains an asymptotic convergence limit. The improved performance 
of LSE for finer grid is expected as the minimization of the residual should satisfy the 
governing equations more accurately. 
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Figure 2.10: Comparison of interpolated and extrapolated pressure with CFD solutions 
for Re = 5.33 with variable lid- velocity. A) Extrapolated pressure compared 
with CFD solution on the extrapolated grid (M 0 ). B) Extrapolated pressure 
compared with CFD solution on fixed fine grid (M„), 121 x 121. (E” : n refers 
to the source of the solution which is compared with the reference solution 
identified by m). 

To follow up on the previous discussion, the variation of the L 2 norm of the 
residual, F, and the L 2 norm error of the extrapolated solution on grid M 0 as compared to 
the fixed grid 121 * 121, with respect to a is examined. Figure 2.1 1 shows that the trend 
of both the accuracy measures are similar but based on minimization of F (Figure 2. 1 1 A) 
LSE picks a value of -1.48 for abut based on L 2 norm error (Figure 2.1 IB) a value of - 
1.6 for a is a better choice. Although it does not identify similar ds that minimize both 
measures, they are close to each other. More assessment will be presented for the flow 
with a constant lid velocity. 

Variable lid-velocity (Re = 1000). The performance of LSE is compared with 
other methods for Re = 1000 and the obtained results are shown in Figure 2.12. The base 
grids used for the extrapolation are 111 x 111 and 121 x 121. It clearly shows in Figure 
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2.12A that for a given fine grid (M 0 ) LSE performs better. Similar to the low Reynolds 
number case, the LSE solution improves as the grid gets finer in terms of the true error 
measure (Figure 2.12B). 


a vs F (Mo= 101) 
ocre 2 = -2.77, log, 0 (F) = -9.945 



A 


a vs E m „(Mo= 101) 
ccre2 = -2.77, log(o(E Mn RE2) = -3.588 



B 


Figure 2. 1 1 : Least square L 2 norm error and error in extrapolated pressure on grid 101 x 
101 for different a for Re = 5.33 and variable lid- velocity. A) Least square L 2 
norm error (F). B) L 2 norm error ( E ” ) (comparing extrapolated pressure with 
CFD solution on grid 121 x 121 (M„)). 

This exercise with variable lid-velocity and two different Reynolds numbers gives 
us incite into few important aspect of LSE. 

• The resolution of the coarse grid has to be adequate for efficient extrapolation. 
Hence for higher Reynolds number, resolution of the base grids is higher. 

• LSE performs better for a considerable range of Reynolds numbers (5.33 and 
1000), thereby showing its efficiency for a wide range of laminar and 
incompressible flow regime. 

• LSE performs better as the grid gets finer since it satisfies the governing equation 
more accurately over the computational domain. 

• The minimization of the residual does not guarantee the minimization of the 
solution error. 
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Figure 2.12: Comparison of interpolated and extrapolated pressure with CFD soluti ons 
for Re = 1000 and variable lid-velocity. A) Extrapolated pressure comp arec j 
with CFD solution on the extrapolated grid (Mo)- B) Extrapolated pressu re 
compared with CFD solution on grid 171 x 171. (£" : n refers to the sou rce 0 f 
the solution which is compared with the reference solution identified by m \ 

Constant lid-velocity (Re = 5.33). At this point it would be interesting to take a 
look at a case where the lid- velocity is uniform to evaluate the efficiency of LSE when 
singularities are present in the solution domain. 

Initially the LSE is implemented on the complete solution domain. In Figure 2.13, 
the sensitivity of L 2 norm of the residual, F and the L 2 norm of the solution error, £* ^ 

with respect to a is shown. LSE identifies a = -0.48 as the optimum value, based on the 
minimization of F (Figure 2. 13 A). But as seen from Figure 2. 13B this value of a is f ar 
from the value required to obtain an accurate solution which is about -4.0. 

To investigate the impact of the singularities, the LSE is carried out over a reduced 
solution domain. Figure 2.14 shows the reduced domain over which the LSE is 
implemented. The extrapolation is carried over a domain with x = 0-> 1 and y = O-^q. 95 . 
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a vs F (Mo = 101) 
ctRE 2 = -2-77, log 10 (F) = -6.282 



aVsE m „ (Mo ■ 101) 

(xre 2 = -2.77, log 10 (E% K ) = -2.365 
aisE = -0.48, log 10 (E M " LSe ) = -2.054 
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Figure 2.13: Least square L 2 norm error and error in extrapolated pressure on grid 101 x 
101 for different a for Re = 5.33 and constant lid-velocity (full domain). A) 
Least square L 2 norm error (F). B) L 2 norm error ( E” ) (comparing 
extrapolated pressure with CFD solution on grid 121 x 121 ( [M n )) . 



Figure 2. 14: The shaded domain is used for LSE. A region between y = 0 0.95 for all * 

is used for extrapolation. 

As shown before in Figure 2.13 for the whole domain, the sensitivity of least 
square error norm, F and the L 2 norm error ( £” ) in comparison to the CFD solution with 
respect to a for the reduced domain, of grid 101 x 101, is shown in Figure 2.15. It is seen 
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that LSE identifies a = -0.48 as the optimum value for a, based on the minimization of F 
(Figure 2. 15 A). But as seen from Figure 2.15B this value of a is not as far from the 
optimum value of a based on the minimization of E ” . The suggested value of a is 

approximately -1.3 as compared to the previous approximate value of -4.0. This suggests 
that a sub-domain which is adequately far from the region of singularity improves the 
performance of LSE. A point to be noted is that during interpolation of data from one 
grid to the other grid, the comer singularity tends to get smeared onto the final grid. Some 
method of treating or accounting for this might have to be looked into. 



A B 


Figure 2.15: Least square L 2 norm error and error in extrapolated pressure on grid 101 x 
101 for different a for Re = 5.33 and constant lid-velocity (reduced domain). 
A) Least square L 2 norm error (F). B) L 2 norm error ( E ” ) (comparing 
extrapolated pressure with CFD solution on grid 121 x 121). 

Figure 2.16 shows the comparisons of different extrapolation techniques. The 
errors shown are based on the reduced domain. There are oscillations noticed in the 
convergence rates as the grid gets finer. This is due to the fact that the region for different 
grids need not always lie between y = 0 -> 0.95, since the nearest grid point may be 
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below or above. This results in different levels of improvement. The present domain 
reduction strategy suggests that when a domain farther away from the singularities is 
considered, the performance of LSE seems better. Of course, this is a very simplistic 
solution and may not necessarily be able to handle other flow problems with more 
complexities. In Figure 2. 16 A, the extrapolated solutions are compared to the fine grid 
(M 0 ) solutions, which show that LSE outperforms the others. In Figure 2.16B, the 
comparison of the extrapolated solution is made with a fixed fine grid solution of size 
121 x 121. The solution on this grid is considered as the benchmark solution. Clearly 
LSE performs well as compared to RE. 


E m n VS N 
Re = 5.33 



Grid Points (N) 

A 
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Figure 2.16: Comparison of interpolated and extrapolated pressure with CFD solutions 
for Re = 5.33 with constant lid- velocity. A) Extrapolated pressure compared 
with CFD solution on the extrapolated grid {M 0 ). B) Extrapolated pressure 
compared with CFD solution on grid 121 x 121 (M n ). (E":n refers to the 
source of the solution which is compared with the reference solution identified 
by m) 

Constant lid-velocity (Re = 1000). LSE is implemented in the same reduced 
domain as shown in Figure 2.14. The coarse grids used are 1 1 1 x ill and 121 x 121 and 
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the fixed fine grid (M„), which is the representative of the correct solution, is 171 x 171. 
Figure 2.17, shows the comparison of different extrapolation schemes. Similar trends are 
noticed for the reduced domain as before. 



Figure 2.17: Comparison of interpolated and extrapolated pressure with CFD solutions 
for Re = 1000 with constant lid- velocity. A) Extrapolated pressure compared 
with CFD solution on the extrapolated grid (M 0 ). B) Extrapolated pressure 
compared with CFD solution on grid 171 x 171. ( E ™ : n refers to the source of 
the solution which is compared with the reference solution identified by m) 

This exercise helps gain an idea of shortcomings related to LSE technique. It shows 
that when singularities are present, the minimization of the L 2 norm of the discretized 
governing equation residual, F, need not minimize the solution accuracy measure, E ” . It 
depends largely on the matrix of equations used for the least square implementation. 
When the singularity is removed and the LSE is carried out on the reduced domain, the 
minimization of F tends towards the minimization of E” and the performance of LSE is 


improved. 
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LSE for coupled pressure-velocity computations 

The coupled algorithm is implement on a lid-driven cavity flow with Re = 5.33 and 
1000 and the variable lid- velocity as given in Eq. (2.36) is used. The geometry is same as 
before (Figure 2.5). The coarse grids {Mi, Mi ) , fine grids (M 0 ) and the fixed grids (M„) 
are all the same as those used for pressure extrapolation. It is noticed that about 3-4 
iterations of this algorithm results in converged values of a, y and with about 3-4 inner 
iterations for the velocity components. Only the error compared to the finest grid (M„) is 
shown since the goal is to obtain a solution as accurate as the exact solution. For clarity, 
only the RE2 and LSE solutions are shown in Figures 2.18 and 2.19 since they performed 
better compared to the rest. 

Re = 5.33. The E” measure for pressure, w-velocity and v-velocity are shown in 

Figures 2. 18 A, 2.18B and 2.18C, respectively. Single mode (constant, ao ) is used for the 
extrapolation of all the flow parameters. The LSE performs well in extrapolating all the 3 
flow parameters onto the finest grid. 

Re = 1000. In Figure 2.19, the LSE solution is compared with the solution of 
second order RE. There is considerable difference in the trends as compared to the low 
Reynolds number case. But as the grid gets finer, the performance of LSE is better than 
RE, suggesting that LSE gives a solution closer to the exact solution than RE. 

This exercise gives the following information: 

• The coupled algorithm is effective in treating all flow variables: pressure, u- and v- 
velocity, adequately, and works well for a wide range of Reynolds number. 

• It is surprising to notice that for the coupled algorithm, the improvement noticed in 
LSE as compared to RE2 seems to be more for the higher Reynolds number case 
t han the lower. Additionally the trend of w-velocity is reversed for Re = 1000 and 
as the grid gets finer, the performance of LSE is poorer in terms of accuracy. This 
requires further study. 
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Figure 2.19: Extrapolated pressure and velocity compared with CFD solution on grid 171 
x 171 ( M „ ), Re = 1000 and variable lid-velocity. A) Pressure. B) w-velocity. 
C) v-velocity. 







CHAPTER 3 

GLOBAL SENSITIVITY AND OPTIMIZATION TECHNIQUES 
In this chapter detail of the surrogate modeling method namely response surface 
methodology (RSM), global sensitivity approach and multi-disciplinary optimization 
(MDO) tools like design of experiments (DOE) and genetic algorithm are presented. 
These tools are applied in the design of a single element injector of a liquid rocket engine 
and the obtained results are presented in chapter 4. 

Response Surface Methodology (RSM) 

The approach of RSM is to perform a series of experiments, based on numerical 
analyses or physical experiments, for a prescribed set of design points (using DOE), and 
to construct a global approximation of the measured quantity (objective) over the design 
space. Figure 3.1 shows the three steps involved for a single-objective optimization study 
which is dependent on two design variables. This is a general case and can be extended 
for multi-objective optimization study as explained later in this chapter. Step I (Figure 
3.1) identifies the designs and the corresponding values of the objective for each design 
over the selected design space. In step II (Figure 3.1) a response surface approximation 
(RSA) of the objective is generated, which is a function of the two design variables. 
Finally step III (Figure 3.1) involves the search for the optimum design using the 
generated RSA for function evaluations. In this dissertation regression analysis is used to 
construct polynomial-based RSAs of assumed order and unknown coefficients. The 
number of coefficients to be evaluated depends on the order of the polynomial and the 
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number of design parameters involved. For instance, a second-order polynomial of N 
design variables has (N+l)(N+2)/(2!) coefficients. A cubic model has 
(N+l)(N+2)(N+3)/(3!) coefficients. In this dissertation, RSAs are constructed using JMP 
(2002), a statistical analysis software that provides a variety of statistical features in an 
interactive format. 

In the practical application of RSM, it is necessary to develop an approximate 
model (RSA) for the true response or the objective. The second order (quadratic) RSA is 
used most frequently as it is the most economical one. Such an approximation for a 
response variable (objective) y with k regressors can be written as 

y = Po+Ya P> x i + Z P» x i + Z Z Pa x i x ) +£ ( 3 -0 

1=1 1*1 i'=l 7=2 

The above equation can be written in matrix notation as follows 
y = Xp + e (3.2) 

where y is the (nx 1) vector of observations, X: (nxn p ) matrix of the levels of the 
independent variables, fi : (n p x 1) vector of the regression coefficients, e: (nx 1) vector of 
random error, n: the number of observations, and n p : the number of terms in the RSA. 

The purpose is to find the vector of least square estimators b that minimizes 

Error = Jf, 2 = e T e = (y - X P) r (y - X p) (3.3) 

1=1 


which yields the least squares estimator of P 
b = (X T X)-'X T y 


(3.4) 
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step I: Define and populate the space 


step II: Generate RSA. 


step III: Search for the optimum 




Figure 3.1: Schematic of the procedure for global design optimization (Papila, 2001). 
The quality of different RSAs can be evaluated by comparing the adjusted rms- 


error a a (Myers and Montgomogery, 1995) value that is defined as: 



where e, is the error at / th point of the data used for fitting. 

The accuracy of the RSA in representing the design space is measured by 
comparing their predictions to the actual values of the objective at test design points, 
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different from those used to generate the fit. The prediction rms-error <7 for the test set is 
given by: 


(3.6) 

In this equation £, is the error at the z th test point and m is the number of test points. 
The test points can be selected using the cross-validation technique (Papila, 2001), which 
is an established method for estimating the prediction accuracy. It is usually performed 
using either a number of random test/train partitions of the data, or a &-fold cross- 
validation (Mullin and Sukthankar, 2000). In &-fold cross-validation, the data points are 
divided in k equally sized subsets. One of the subset is kept for testing and the remaining 
k - 1 are used to fit the RSA. This is done for all the k subsets and the average error of all 
the k RSAs is computed. The fitting £-1 sets which has the minimum prediction rms- 
error, a, is selected and the remaining subset used as the testing set. 

In certain cases, especially when higher order polynomials are used, the data 
available for RSA may not be enough to spare some for testing the fitted model. Hence, 
an alternate method to estimate the performance of the RSA is to compute the PRESS 
rms-error. This method was proposed by Allen (1971, 1974) and it computes a sum of 
squares of the residuals. The residual is obtained by fitting a RSA over the design space 
after dropping one design point from the fitting set and then comparing the predicted 
value of the RSA for that point with the expected value. This is done for every point in 
the design space. The PRESS rms-error is given by 



PRESS 


1 


Eb .-£] 2 


/=! 


n 


(3.7) 
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where y,- is the expected value, y ( . is the value predicted by the RSA for the I th point, 
which is excluded while generating the RSA and n is the number of design points. If this 
value is close to the adjusted rms-eror, cr a , then the model performs well. 

A measure of the variability in an objective accounted by its RSA is given by the 
coefficient of multiple determination (Myers and Montgomogery, 1995), R , given as 

r2=1 _^e_ ( 3 . 8 ) 

SS 

yy 

n 

where SSe is the sum of squares of the residuals (=^(y,. -y,.) 2 ) where y is the predicted 

»=i 

value by the RSA. SSyy is the total sum of squares about the mean given by 

ss„=ss s +ss R =Y,(y,-y? 

i=J 

where y is an overall average ofy,. SSr is the sum of squares due to regression 


n 

(= ^ ( j>. -y ) 2 ) where y is the overall average of y,-. Since R 2 increases as more terms are 

i=i 

added in the RSA, the adjusted coefficient of multiple determination R 2 , a normalized 
measure is usually preferred. It accounts for the degrees of freedom in the model and is 


given by 



SS E /(n-n p ) 
SS yy /(n- 1) 


= 1 - 



(3.10) 


For a good fit, the value of R] should be closer to one. 

The polynomial-based RSA is effective in representing the global characteristics of 
the design space. It can filter the noise associated with design data. The linearity of the 
polynomial-based RSA allows us to use statistical techniques known as design of 
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experiments (DOE) to find efficient fitting sets. On the other hand, depending on the 
order of polynomial employed and the shape of the actual response surface, the RSA can 
introduce a substantial error in certain region of the design space. 

Design of Experiments (DOE) 

It is well established that the predictive capability of RSM is greatly influenced by 
the distribution of sampling points in design space (Unal et al., 1997, 1998). Design of 
experiment is required to select designs for RSA that minimizes the effect of noise. There 
are different types of DOE techniques in the literature as reported by Haftka et al. (1998). 
One of the most conservative DOE techniques available in literature is the full factorial 
design (JMP, 2002). Unal et al. (1996, 1997) studied response surface modeling using 
orthogonal arrays (O A) in computer experiments for reusable launch vehicles and 
illustrated that using this technique can minimize design, development, test and 
evaluation cost. Unal and Dean (1995) studied a robust design method based on the 
Taguchi method (Unal and Dean, 1991; Dean, 1996) to determine the optimum 
configuration of design parameters for performance, quality and cost. They demonstrated 
that using such a robust design method for selection of design points is a systematic and 
efficient approach for determining the optimum configuration. The full factorial design 
and the OA are explained below. 

In full factorial design the range of each design variable is divided into one or more 
intervals, which mark the number of levels. These intervals are usually, evenly spaced. 
All the possible combinations of the levels of all the design variables give the design 
points in the design space. For example, for three variables with three levels each there 
are totally 27 design points (Figure 3.2). 



50 



Figure 3.2: Full factorial 3-level design (dark circles refer to points in the foreground and 
light circles refer to points in the background). 

An OA is a fractional factorial matrix that assures a balanced comparison of levels 
of any factor or interaction of factors. Consider A, a matrix with elements of aj where j 

denotes the row (j = 1,2... n r ) and i denotes the column (i = l,2...n c ) that aj belongs to, 
supposing that each aj e Q = {0, 1 . . .q- 1 } . If the columns representing each design 
variable are mutually orthogonal, matrix A is called an orthogonal array. It is of strength t 
<n c if in each n r -row-by-t-column sub-matrix, the q‘ possible distinct rows occur X (index 
of the array) times. Such an array is denoted by OA(n r ,n a q,t ) by Owen (1992). Since the 
points are not necessarily at the vertices, the OA can be more robust than the full factorial 
design in interior design space and is less likely to fail the analysis tool. Based on the 
DOE theory, OA can significantly reduce the number of experimental configurations 
(Papila, 2001). In this dissertation OA is used to obtain the design for the optimization 
study. 
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Sensitivity Analysis 

Unlike local sensitivity where the partial derivatives are used to locally estimate the 
sensitivity of an objective to a specific design variable, global sensitivity allows the study 
of overall model behavior. To understand the concept, assume a surrogate model (in this 
dissertation a RSA) of an objective, J{X) as a function of the design variables, X, where X 
is the vector of the design variables xi, X 2 , ..., x„, scaled between zero and one. This 
surrogate model (in this dissertation, polynomials of n* order) represented in Eq. (3.11) 
has to be square integrable. 


f{,X ) — fo + ^ fi (*, ) + ^ fij (-*, , JCy )+ ... + f\ 2 ...n (*1 ’ X 2 >"’> X n ) 


(3.11) 


«</ 


An analysis of variance (ANOVA) (Archer et al., 1997) representation of Eq. (3.1 1) 
can be obtained by imposing the following condition 

=0 (3.12) 

0 


for k = ii , ..., i s . This condition is essential for the uniqueness of Eq. (3.11). This results 
in a set of summands of different dimensions. Each of the components is responsible 


for the interactive contribution of design variables ,...,x fj to the variability of the 

objective flX) over the n-dimensional unit hypercube design space. 

Integrating Eq. (3.1 1) over all the design variables gives 
j/(x)n<& =/„ (3.13) 

Integrating Eq. (3.1 1) over all the design variables except x, gives 

f/«n d ** = f & + ^ ( x i ) < 3 - 14 ) 

k*i 


from which/(x,) is obtained. 
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Integrating Eq. (3.11) over all the design variables except jc, and xj gives 
|/Wn^=/o +fiM+fj(xj)+f,j(x lt x j) (3.15) 

k*i,j 

from which fjixuxj ) is obtained. This can similarly be extended to obtain the remaining 
summands offlX). The imposed condition ensures that the summands are orthogonal in 
nature and since/(X) is square integrable the summands are too. Therefore the partial 
variances can be calculated as 


D, i = f fj 2 1 dx ...dx i 

l i— l s J*' 1 !— 1 ! *1 l t 


and the total variances is 


(3.16) 


D= jf 2 (x)dx-f 0 2 (3.17) 

Therefore it can be shown (Sobol, 1993) that 


n n 

0 = L 2X..,. 


ij <...</, 


(3.18) 


Partial variance gives a measure of the spread of the function due to one of the 
independent variables. Total variance gives a measure of the spread of the dependent data 
due to all the independent variables. Sensitivity is a measure of the contribution of an 
independent variable to the total variance of the dependent data. Sobol (1993) had 
proposed a variance-based non-parametric approach to estimate the global sensitivity 
using Monte Carlo methods. To calculate the total sensitivity of any design variable, say 
xi, the design variable vector X is divided into two complementary subsets, xi and Z 
where Z is a vector containing^, X 3 , ..., x„. The purpose of using these subsets is to 
isolate the influence of x, on f(X) variability from the influence of the remaining design 
variables included in Z. The total sensitivity index for*/ is then defined as 
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(3.19) 

where 

D? = Z>,+Z>,„z (3.20) 



D x is the partial variance associated with x; and D x Z is the sum of all partial variances 

associated with any combination of the remaining variables representing the interactions 
between xj and Z. Similarly the partial variance associated with Z can be defined as Dz- 
Therefore the total objective variability can be written as 

D-D^D Z +D^ (3.21) 

In cases where the/(X) is an analytical function, the multidimensional integrals for 
computing the partial variances can be evaluated analytically. However, the Monte Carlo 
approach is applicable under general conditions (e.g., any model, design under 
uncertainty) and has been adopted in this dissertation. The designs for the Monte Carlo 
approach are selected without any preference of one design over the other from the 
uniformly distributed design space of unit sides. Hence the design variables are 
uncorrelated which is essential for the effective implementation of the method proposed 
by Sobol (1993). The variance estimates can be obtained using the following expressions: 

(322) 

where fo is the mean of the objective. 


* ’ y=l 


(3.23) 


JV j ={ 


(3.24) 



where the terms (xj, Z) and (xi ’, Z’) represent random sample designs. Equations. (3.23), 
(3.24) and (3.25) give the estimates for D Xi , D and Dz, respectively and D z can be 

calculated using Eq. (3.21). Once these estimates are known the total sensitivity index of 
the objective variability with respect to a given design variable can be obtained using Eq. 
(3.19). The influence of a design variable to an objective variability without accounting 
for any of its interactions with other variables is denoted as a main factor index and given 
as 

S„=DJD (3.26) 

Each pair of random samples requires three different objective function evaluations 
(e.g .,J{x], Z ),/ (xj, Z’) and / (xj Z)). The mean and the total variance of an objective 
need to be estimated only once, and only two Monte Carlo integrals per design variable 
are necessary to compute the main factor and total sensitivity indices. Hence the increase 
in computational cost is linear with the increase in the number of design variables. 

These indices can be used to understand the influence of design variables on the 
variability of any given objective over the chosen design space. This method proposed by 
Sobol (1993) is effective when the design variables are uncorrelated. In a multi-objective 
optimization study multiple optima are obtained which collectively give the Pareto 
optimal front (POF). The design variables at different regions along the POF share some 
similar features. Therefore, the design variables are correlated and this correlation has to 
be accounted for before estimating the influence of a design variable on an objective. The 
total sensitivity and main factor indices, which require uncorrelated data, cannot be used 
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for such a situation. Hence, for this purpose, a partial correlation coefficient (JMP, 2002) 
was calculated. 

Estimation of partial correlation coefficient (JMP, 2002) is a variance-based 
parametric approach that gives a measure of the linear relationship between the variances 
of a design variable, say xi and an objective J{X), after the influence of other variables 
have been filtered out. It gives a measure of expected change in_/(A) per unit change in xj 
which in other words gives the sensitivity of the function to the design variable. Linear 
approximations are obtained for*; and/(X) as a function of the components of Z and the 
residuals (say, ri and r 2 , respectively) measured as the difference between the data used 
for the approximations and the predicted value of the approximations are estimated. A 
partial correlation coefficient is the correlation between these two residuals and is given 
by 

fKvJiKvJi)] (3 . 27) 

<J cr 

r l r 2 

Optimization Algorithm 

The optimization process can be divided into two parts: 

1. RS A phase, 

2. Optimizer phase. 

In the first phase, RSA are generated with the available data set. In the second 
phase the optimizer uses the RSA during the gradient-based search for the optimum (the 
maximum or minimum of the objective). The initial design for the search is randomly 
selected from within the design space. The flowchart of the process is shown in Figure 
3.3. The optimization problem at hand can be formulated as min {/(X)} subject to lb < X < 
ub, where lb is the lower boundary vector and ub is the upper boundary vector of the 
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design variables vector X. If the goal is to maximize the objective function then/fx) can 
be written as -g(X), where g(X) is the objective function. Additional linear or nonlinear 
constraints can be incorporated if required. Solver (2002), an optimization tool available 
as part of Microsoft Excel package is an ideal tool for such studies. This tool uses the 
Generalized Reduced Gradient ( GRG2 ) nonlinear optimization code developed by 
Lasdon et al. (1978). 


r j 



Figure 3.3: Two phases of the optimization process, where Phase 1 deals with RSA and 
Phase 2 deals with optimization. 

In most optimization studies it is desirable to simultaneously optimize more than 
one response/objective. One method is to build, from the individual objectives, a 
composite objective known as the desirability function. The method allows for the 
designers’ priorities to be addressed during the optimization procedure. The first step is to 
define the desirability function, d for each objective. Desirability function is a 
normalizing of the objective such that its’ maximum and minimum over the design space 
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lie between zero and one. In the case where an objective should be maximized, say Ri, 
the desirability takes the form: 

d.= (3.28) 

\B-A) 

where B is the target value and A is the lowest acceptable value such that dj = 1 for any 
Rj> B and di = 0 for Rj < A. The weighting factor t is set according to the designer’s 
goal and the compromise he wants to achieve between different objectives. 

In the case where a response is to be minimized, say R 2 , the desirability takes on 
the form: 



(3.29) 


where C is the target value and E is the highest acceptable value such that d 2 = 1 for any 
R 2 <C and d 2 = 0 for R 2 > E. Choices for A, B, C, and E are made according to the 
designer’s priorities or simply as the maximum and minimum values of the objective 
over the design space. Values of s and t are set based on the priority of the objective. The 
sensitivity of the parameter, s, illustrated in Figure 3.4 can be instructive (Myers and 
Montgomery, 1995). Desirability functions with s « 1 imply that a product need not be 
close to the response target value, C, to be quite acceptable. However, s = 8, implies that 
the product is nearly unacceptable unless the response is close to C. A single composite 
desirability D is defined based on the geometric mean of the desirabilities of the 
individual objectives. 

D = {d r d 2 -d 3 ....d m f"> (3.30) 

This is then submitted to an optimization toolbox to be maximized. 
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Figure 3.4: Desirability function (d) for various weight factors, s. (Note: B < A) (Myers 
and Montgomery, 1995). 

Multi-Objective Genetic Algorithm 

One of the other optimization approaches is to use Multi-objective genetic 
algorithm (MOGA). It works on the principle of survival of the fittest. Genetic operators 
like reproduction, crossover and mutation are employed to find the optimal solution. In a 
multi-objective optimization scenario, when the objectives are conflicting in nature, many 
representative optimal solutions can be obtained. All these solutions comprise a set called 
the Pareto optimal set (Miettinen, 1999). The solutions in this set are considered non- 
dominated as they are not inferior to any other solution in the entire design space without 
having a bias towards some of the objectives. The function space of all the solutions in 
the Pareto optimal set is termed the Pareto optimal front (POF) (Miettinen, 1999). When 
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there are two objectives, the POF is a curve. When there are three objectives, the POF is 
represented by a surface. If there are more than three objectives, it is represented by 
hyper-surfaces. 

One of the MOGA that has been shown to be effective in finding the Pareto- 
optimal solutions is the elitist non-dominated sorting genetic algorithm (NSGA-II) (Deb 
et al., 2000). The algorithm can be described as: 

1 . Randomly initialize population (designs in the design space) of size N. 

2. Compute objectives for each design. 

3 . Rank the population using non-domination criteria. 

4. Compute crowding distance (this distance finds the relative closeness of the 
solution with other solutions in the objective space.) 

5. Employ genetic operators - selection, crossover & mutation - create new 
population. 

6. Evaluate objectives for the new population. 

7. Combine the two populations, rank them and compute the crowding distance. 

8. Select N best individuals. 

9. Go to step 3 and repeat till termination criteria is reached, which in the 
current study is chosen to be the number of generations. 

It is shown in a number of studies that using a combination of MOGA and local 
search (also known as hybrid GA), helps achieve faster convergence to the global Pareto 
optimal solution set (Deb and Goel, 2000, 2001; Goel, 2001; Goel and Deb, 2002; 
Ishibuchi and Yoshida, 2002). The a posteriori hybrid method (Deb and Goel, 2001) 
assumes the set of solutions generated by MOGA simulations as a starting point for the 
local search. Most of the local search methods are very efficient only for single objective 
optimization problems. Hence, a e-constraint strategy (Deb and Goel, 2001) helps reduce 
the dimensionality of the problem. Sequential quadratic programming, available in 
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Matlab (2002) can be used as the local search optimizer. This gives a set of solutions 
from which the dominated and duplicate solutions are removed to obtain the global 
Pareto optimal solution set and the POF. 

Hierarchical Clustering Method 

The POF can be divided into a number of clusters using a hierarchical clustering 
algorithm (Jain and Dubes, 1988) to assist the designer in selecting the optimal solution 
of choice. The clustering algorithm can be summarized as: 

1 . Start by assuming all the solutions as individual clusters. 

2. Find the mean of each cluster. 

3. Find the distance between clusters. 

4. Merge closest clusters. 

5. Go to step 2 till the number of clusters is equal to some prefixed value. 

6. Find the member of each cluster closest to the mean of the cluster. This is the 
representative element. 

Visualization Using Box Plot 

Box plot (JMP, 2002) is a visualization tool, which can help understand the 
variability in the design variables and objectives within the clusters of the POF. It can 
also assist in identifying the design variable that should be fixed when analyzing a given 
cluster. A schematic of a box plot with the y-axis representing the range of the scaled 
design variable is shown in Figure 3.5. 

Figure 3.5 shows that 25%, 50% and 75% of the data lie below the lower hinge, 
median and upper hinge, respectively. The difference between the upper and lower hinges 
is known as the “H-spread”. The inner fence is located 1 step beyond the hinges, which is 
equal to 1.5 times the H-spread. The upper adjacent value is identified as the largest value 
below the upper inner fence. The lower inner fence and lower adjacent value can be 
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similarly determined. These box plots will be used to visualize the variability along the 
POF. The spread between the upper and lower adjacent values gives the range of 
variation of a variable in a given set. Any variable lying beyond this spread is a potential 
outlier of the set. Small range gives a tight bound on that particular variable. Sometimes 
the box plot is known to collapse to a single point, which suggests that the variable 
should be fixed at that value. 



Figure 3.5: Schematic of a box plot 




CHAPTER 4 

DESIGN OPTIMIZATION - SINGLE ELEMENT ROCKET INJECTOR 

In this chapter the sensitivity and optimization studies of a single element liquid 
rocket injector design are presented. Initially a literature survey of the past injector design 
studies is presented. Following this the injector model is described. Then the results 
related to the sensitivity analyses and optimization studies are presented. 

Literature Review 

A critical goal for space propulsion design is to make the device safer, more 
affordable and more reliable. The design of combustion devices, namely, injectors, 
chambers and nozzles, will help in meeting these goals. The characteristics of the injector 
design are a key factor for both performance and thrust chamber environments. The thrust 
chamber performance is estimated by the rate and the extent to which mixing and 
resultant combustion occurs. The location of the mixing and resultant combustion 
determines the injector and thrust chamber thermal environments. These environments 
include temperatures on the combustor wall, the injector face and, for coaxial injectors, 
the oxidizer post tip. The difficulty encountered in designing injectors that perform well 
and have manageable environments is that the factors that promote performance often 
lead to increased heating of the solid surfaces of the injector and combustor thereby 
reducing the life expectancy or survivability of these components. 

Current injector design tools have been in use for 30 years or more and are largely 
empirical-based (Rupe, 1965; Dickerson et al., 1968; Pavli, 1979; Nurick, 1990; Pieper, 
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1991). Extensive sub- and full-scale hot-fire test programs often guided these injector 
design methodologies. The experimental databases, and thus the tools developed from 
them, are limited, in terms of design space, to specific element configurations that have 
been tested (Calhoon et al., 1973). In terms of scope, the design tools typically focus on 
performance, with the environment being a secondary consideration. The limited amount 
of environmental information available from these tools is usually one-dimensional and 
not functionally related to details of the injector design. It is very doubtful that 
application of these traditional design tools will result in robust future propulsion devices. 
Over time space propulsion programs with compressed schedules, lower budgets and 
more stringent requirements have resulted in the development of broader and more 
efficient injector design methodologies. 

The initial work by Tucker et al. (1998) and Vaidyanathan et al. (2000) focused on 
the optimization of a shear coaxial injector element with gaseous oxygen (GO 2 ) and 
gaseous hydrogen (GH 2 ) propellants. In this study, the design data was generated using 
an empirical design methodology developed by Calhoon et al. (1973). Calhoon et al. 
(1973) conducted a large number of cold-flow and hot-fire tests over a range of 
propellant mixture ratios, propellant velocity ratios and chamber pressure for shear 
coaxial, swirl coaxial, impinging, and premixed elements. The data were correlated 
directly with injector/chamber design parameters, which were recognized from both 
theoretical and empirical standpoints as the controlling variables. For the shear coaxial 
element, performance, as measured by energy release efficiency, ERE, was obtained 
using correlations taking into account combustor length, L CO mb (length from injector to 
throat) and the propellant velocity ratio, V/V a . The nominal chamber wall heat flux at a 
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point just downstream of the injector, Q„ 0 m, was calculated using a modified Bartz 
equation and was correlated with propellant mixture ratio, O/F, and propellant velocity 
ratio, V/V 0 to yield the actual chamber wall heat flux, Q. The objective was to maximize 
injector performance while minimizing chamber wall heat flux (lower heat fluxes reduce 
cooling requirements and increase chamber life/survivability) and chamber length 
(shorter chambers lower engine weight). 

In Tucker et al. (2000) the designs of an impinging injector element and a swirl co- 
axial injector element have been carried out. For the impinging injector element, the 
empirical design methodology of Calhoon et al. (1973) used the oxidizer pressure drop, 
AP 0 , fuel pressure drop, APf, combustor length, L com b, and the impingement half-angle, a 
as design variables. Objectives included ERE (a measure of element performance), wall 
heat flux, Q w , injector heat flux, Q inJ , relative combustor weight, W re i, and relative injector 
cost, C rd . The gaseous propellants were injected at a temperature of 540R. The empirical 
design methodology used to characterize the ERE and Q w used a quantity called the 
normalized injection momentum ratio to correlate the mixing at the different design 
points for the triplet element. 

The swirl coaxial element has been used sparingly in USA, but has been widely 
used in Russia because of its reported ability to perform well over a large throttle range 
(Gill and Nurick, 1976). This has sparked the interest in exploring its possibilities. The 
empirical design methodology of Calhoon et al. (1973) used the oxidizer pressure drop, 
AP 0 , fuel pressure drop, APf, , combustor length, L com b, and the full cone swirl angle, 0, as 
design variables. The objectives modeled were ERE (a measure of element performance), 
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wall heat flux, Q w , injector heat flux, Q in j, relative combustor weight, W re i, and relative 
injector cost, C re i- 

The advent and advancement of CFD in the last 20 years have produced a 
capability that has shown potential as an improved design tool in many areas of rocket 
propulsion. The application of CFD to injector design has lagged behind other areas such 
as turbomachinery because the physical models are more complicated for multiphase, 
turbulent reacting flows. New models that efficiently account for some of the complex 
processes (Cheng and Farmer, 2002) and thus increase the solution fidelity have recently 
become available. However, the three-dimensional geometry of multi-element injectors 
and the complex physical processes inherent in the flows that issue from them create 
major obstacles in validating the solution and generating significant amount of parametric 
solutions. The harsh high pressure and temperature environments typical of injector flows 
create significant difficulty in obtaining experimental data of satisfactory quality to 
validate and guide further development of computational models. Finally, solving the 
equations, for multiphase reacting flows, with high resolution typically requires lengthy 
computational times. However, continuing increases in computer speed and progress in 
parallel processing have begun to mitigate this turnaround problem. 

It has long been known that small changes in injector geometry can have significant 
impact on performance, as well as on environmental variables such as combustion 
chamber wall and injector face temperatures and heat fluxes (Gill and Nurick, 1976). 
Several geometric design variables can be accounted over appropriate ranges by 
calculating the performance and environmental variables of the injector using CFD. 
Global approximation method like RSM can provide a continuous representation of the 
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objectives over the design space thereby reducing the number of numerical computations 
required for sensitivity analyses. This global approximation can also guide the 
optimization study. Additionally, a global approximation method can also at times 
address issues related to CFD analysis, for e.g. grid resolution, thereby increasing the 
fidelity of the computations. The present work is first of its kind where a CFD-based 
design optimization methodology is proposed for the design of an injector (Vaidyanathan 
et al., 2003a; Vaidyanathan et al., 2004b). 

Scope 

The focus is on sensitivity and trade-off analyses for the design of a gaseous 
injector for liquid rocket propulsion. The data for the analyses is obtained from surrogate 
models (RSAs) of the objectives, and the multi-objective optima (Pareto optimal front, 
POF) generated by Goel et al. (2004) with the aid of multi-objective genetic algorithm 
(MOGA) and a local search method. The regions of the POF that represent different 
trade-offs among the objectives are obtained through a hierarchical clustering algorithm. 
The initial study concentrates on the generation of response surface approximations 
(RSAs) and preliminary optimization studies. Then it is followed by an elaborate 
sensitivity and trade-off analyses. The later part of the study is broadly divided into two 
parts. Firstly a sensitivity study is carried over the whole design space. The contribution 
of the design variables to the objective variability is calculated using a variance-based 
non-parametric approach (Sobol, 1993) and correlations between objectives are 
investigated. Secondly, sensitivity analyses are conducted on clusters along the POF. Box 
plots are used to highlight the variability of the design variables and the objectives within 
a given cluster. Additionally, the linear relationships between the variance of the design 
variables and the objectives are explored with the aid of partial correlation coefficients. 
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Injector model 

Liquid rocket propulsion injector elements can be categorized into two basic types 
based on propellant mixing. The first type is an impinging element (Figure 4.1 A) where 
mixing occurs by direct impingement of the propellant streams at an acute angle. The 
impingement enhances mixing by head-on interaction between the oxidizer and fuel (Gill 
and Nurick, 1976). The second type of injector consists of non-impinging elements where 
the propellant streams flow in parallel, typically in coaxial fashion (Figure 4. IB). Here, 
mixing is accomplished through a shear-mixing process (Calhoon et al., 1973). From a 
design standpoint, both element types have appealing characteristics. However both also 
have undesirable characteristics. For instance, if the impinging element has an F-O-F 
arrangement (Figure 4.1 A), the mixing occurs rapidly, which can yield high performance. 
However, since the combustion zone is close to the injector face, the potential for high 
levels of injector face heating must be considered. If the non-impinging element is 
assumed to be a shear coaxial element, mixing across the shear layer is relatively slow, 
requiring longer chambers to allow complete combustion. However, since the combustion 
zone is spread over a longer axial distance, the injector face is generally exposed to less 
severe thermal environments. 

Important design parameters for the impinging element (assuming fixed mass flow 
rates and constant propellant inlet conditions) include relative orifice size (or, relative 
stream momentum ratio), impingement angle and orifice spacing. Important parameters 
for the shear coaxial element (assuming fixed mass flow rates and constant propellant 
inlet conditions) include the area ratio of the two concentric tubes (or velocity ratio) and 
the shear area between the two propellant streams (i.e., the oxidizer post tip thickness) 
(Calhoon et al., 1973). 
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A B 


Figure 4.1: Schematic of the G02/GH2 impinging and coaxial injector elements. A) F-O- 
F impinging element. B) Coaxial element. 

One can combine the above-mentioned characteristics of injector types to develop 
hybrid concepts. For example, it has been noted (Gill and Nurick, 1976) that performance 
improvement in the shear coaxial element can be realized by directing the fuel toward the 
oxidizer stream rather than parallel to it and thinning the oxidizer post wall. The first 
modification causes the shear coaxial element to take on some of the aspects of the F-O-F 
impinging element. These notions lead to the hybrid element shown in Figure 4.2 that has 
been developed by The Boeing Company (U. S. Patent 6253539). 



Figure 4.2: Schematic of Hybrid Boeing Element (U. S. Patent 6253539) 


The four independent design variables chosen for the element used in this study are 


shown in Figure 4.3 A. They are the angle at which the H 2 is directed toward the oxidizer, 
a , the change in H 2 flow area from a baseline area, AHA, the change in 0 2 flow area from 
a baseline area, AO A, and the oxidizer post tip thickness, OPTT. The fuel and oxidizer 



















flow rates are held constant. The independent variable ranges considered in this exercise 


are shown in Table 4. 1 


0 2 Post Tip Thickness 
(OPTT=x-2x in) 


O 2 Flow Area 
(AOA=0~40%) 


H 2 Flow Area 
(AHA =0-25%) 


Figure 4.3: Design variables and objectives of the single element rocket injector. A) 
Design variables and their range. B) Objective functions. 
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Table 4.1: Range of design variables (a is an acute angle in degrees and x is the thickness 


of OPTT in inches) 

Variable 

Minimum 


Maximum 



Actual 

Scaled 

Actual 

Scaled 

a 

0° 

0 

a 0 

1 

AHA 

Baseline 

0 

Baseline+25% 

1 

AOA 

Baseline-40% 

0 

Baseline 

1 

OPTT 

x in 

0 

2x in 

1 


The dependent variables chosen for design evaluation are shown in Figure 4.3B. 
They are the maximum injector face temperature, TF max , the wall temperature at a 
distance three inches from the injector face, TW 4 , the maximum oxidizer post tip 
temperature, TT max and centerline axial location where the combustion is 99% complete, 
(hereafter referred to as the combustion length) Xc c (Figure 4.3B). The three temperatures 
(calculated as adiabatic wall temperatures in this study) were chosen as indicators of local 
environments. Lower temperatures would indicate a design that had longer 
life/survivability due to decreased thermal strain on the part. The combustion length, Xc C , 
was chosen as a measure of performance. Shorter combustion lengths indicate better 
performance. 

Numerical Procedure 

A pressure-based, finite difference, Navier-Stokes solver, FDNS500-CVS (Chen, 
1989; Wang and Chen, 1990; Chen and Farmer, 1991) is used in this study. The Navier- 
Stokes equations, the two-equation turbulence model, and kinetic equations are solved. 
Convection terms are discretized using either a second order upwind, third order upwind 
or a central difference scheme, with adaptively added second order and fourth order 
dissipation terms. For the viscous and source terms, second order central difference 
scheme is used. First order upwind scheme is used for scalar quantities, like turbulence 
kinetic energy and species mass fractions, to ensure positive values. Steady state is 
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assumed and an implicit Euler time marching scheme is used for computational 
efficiency. The chemical species continuity equations represent the H 2 - O 2 chemistry. It 
is represented with the aid of a 7-species and 9-reaction set (Chen, 1989; Wang and 
Chen, 1990; Chen and Farmer, 1991). The simulation domain and the boundary 
conditions used in all the CFD cases are shown in Figure 4.4. Because of the very large 
aspect ratio, both the injector and chamber have been shortened (at the cross hatched 
areas) for clarity. Both fuel and oxidizer flow in through the west boundary where the 
mass flow rate is fixed for both streams. The nozzle exit, at the east boundary, is modeled 
by an outlet boundary condition. The south boundary is modeled with the symmetry 
condition. All walls (both sides of the oxygen post, the outside of the fuel annulus, the 
outside chamber wall, and the faceplate) are modeled with the no slip adiabatic wall 
boundary condition. Each CFD analysis was done on a 200 CPU PC cluster with an 
AMD Athlon MP 1800 (1.8 GHz) chip and 1 GB RAM. All the cases were run 
concurrently and took about 5 days. 


adiabatic walk 


injector face 


oxygen po st 



'symmetry 

Figure 4.4: Simulation domain and boundary conditions 

Design Space 

In this study OA (with n c = 4, q = 2,t = 3 and X = 2) is used to generate 54 designs 


for RSA and testing of the approximation. To test the RSA, 14 of the 54 designs are 


72 


selected using cross-validation techniques (Papila, 2001). During the CFD simulations, 
two of the 40 designs (fitting set) were found to be unacceptable because they exhibited 
unsteady behaviors while the numerical algorithm used was a steady state model. Hence, 
the final information included 38 designs for fitting the RSA and 14 to test their 
predictive capabilities. All the design variables are scaled between zero and one based on 
their upper and lower bounds. All the objectives obtained from the CFD solutions of the 
52 valid designs are scaled between zero and one based on the maximum and minimum 
values and polynomial-based RS As generated. Once the RSAs are generated, the scaled 
objectives are normalized between zero and one based on the maximum and minimum of 
the generated RSAs, which will be referred to as the normalized objectives in the text. 
The scaled values are the non-normalized values of the objectives. The fitting and testing 
designs are shown in Tables 4.2 and 4.3, respectively. It should be noted that when the 
AOA is 1 or 0, the O 2 flow area is reduced by 0% or 40%, respectively, as compared to 
the baseline area. 

Results and Discussion 

CFD Analysis 

Comparison of two of the evaluated designs serves to illustrate the motivation for 
combining CFD analyses and an efficient global approximation model during the design 
process. The scaled independent design variables are shown for the two cases in Table 
4.4. In terms of the design space evaluated, these two designs are seen to be quite 
different. The normalized dependent variables are also shown in Table 4.4 with the 
temperatures shown in contour plots (Figures 4.5 and 4.6). 
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Table 4.2: Fitting designs (unacceptable designs are marked in bold). 


Case # 

a 

AHA 

AOA 

OPTT 

Case 1 : 

0 

0 

1 

0 

Case 2: 

0 

0 

1 

0.5 

Case 3: 

0 

0 

1 

1 

Case 4: 

0.5 

0.5 

0 

1 

Case 5: 

1 

1 

0.5 

0 

Case 6: 

1 

1 

0.5 

0.5 

Case 7: 

1 

1 

0.5 

1 

Case 8: 

0 

0.5 

0.5 

0 

Case 9: 

0.5 

1 

1 

0 

Case 10: 

0.5 

1 

1 

0.5 

Case 11: 

0.5 

1 

1 

1 

Case 12: 

1 

0 

0 

0 

Case 13: 

1 

0 

0 

0.5 

Case 14: 

1 

0 

0 

1 

Case 15: 

0 

1 

0 

0 

Case 16: 

0 

1 

0 

0.5 

Case 17: 

0 

1 

1 

1 

Case 18: 

0.5 

0 

0.5 

0 

Case 19: 

0.5 

0 

0.5 

1 

Case 20: 

1 

0.5 

1 

0 

Case 21: 

1 

0.5 

1 

0.5 

Case 22: 

1 

0.5 

1 

1 

Case 23: 

0 

1 

0.5 

0 

Case 24: 

0 

1 

0.5 

1 

Case 25: 

0.5 

0 

1 

0 

Case 26: 

0.5 

0 

1 

1 

Case 27: 

1 

0.5 

0 

0 

Case 28: 

1 

0.5 

0 

1 

Case 29: 

0 

0 

0 

0 

Case 30: 

0 

0 

0 

0.5 

Case 31: 

0 

0 

1 

1 

Case 32: 

0.5 

0.5 

0.5 

0.5 

Case 33: 

1 

1 

1 

0 

Case 34: 

1 

1 

1 

1 

Case 35: 

0 

0.5 

1 

0 

Case 36: 

0 

0.5 

1 

1 

Case 37: 

0.5 

1 

0 

0 

Case 38: 

0.5 

1 

0 

1 

Case 39: 

1 

0 

0.5 

0 

Case 40: 

1 

0 

0.5 

1 
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Table 4.3: Testing designs.. 


Case # 

a 


AOA 

OPTT 

Case 41: 

0.5 

0.5 

0 

0 

Case 42: 

0.5 

0.5 

0 

0.5 

Case 43: 

0 

0.5 

0.5 

0.5 

Case 44: 

0 

0.5 

0.5 

1 

Case 45: 

0.5 

0 

0.5 

0.5 

Case 46: 

0 

1 

0.5 

0.5 

Case 47: 

0.5 

0 

1 

0.5 

Case 48: 

1 

0.5 

0 

0.5 

Case 49: 

0.5 

0.5 

0.5 

0 

Case 50: 

0.5 

0.5 

0.5 

1 

Case 5 1 : 

1 

1 

1 

0.5 

Case 52: 

0 

0.5 

1 

0.5 

Case 53: 

0.5 

1 

0 

0.5 

Case 54: 

1 

0 

0.5 

0.5 


Table 4.4: Independent and dependent variable (objectives) for Cases X and Y from CFD 


computations (non-normalized values shown). 


Case 

a 

AHA 

AOA 

OPTT 

TF ma x 

tw 4 

TT m ax 

Xcc 

X 

1 

0 

0 

0 

0.998 

0.927 

0.128 

-0.004 

Y 

0 

0.5 

0.5 

1 

0.285 

0.395 

0.923 

0.567 



Figure 4.5: Temperature field for Cases X and Y. 


The chamber wall and injector face temperatures for Case Y (as seen in Figure 4.5) 
are low or moderately low, while for Case X, they are high. Figure 4.7 shows a large 
recirculation zone located between the injector and the chamber wall. This recirculating 
flow strips hot gases from the flame and causes them to flow back along the chamber 
wall and injector face. This phenomenon regulates the chamber wall temperature, TW4 
and the injector face temperature, TF max . Figure 4.6 shows that the other life/survivability 
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indicating variable, the maximum oxidizer post tip temperature, TT max , has essentially the 
opposite trend as compared to the other two temperatures. 



Figure 4.6: Near-injector temperature field for Cases X and Y. 



Large 

recirculation zone 


0 2 post tip 


wall 


Injector face 


Figure 4.7: Large recirculation zone in the combustion chamber. 
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The performance indicator, combustion length, X,* is seen (Table 4.4) to be at the 
minimum for Case X (shorter combustion lengths indicate better mixing elements) and at 
a moderate level for Case Y. Given these observations, it is clear that the dependent 
variables exhibit competing trends such that no design will produce the “best” values for 
all the dependent variables as desired in this study. 

These comparisons confirm the past observation that changes in the injector design 
details have major effects on injector performance and injector-generated environments. 
Injector designs which address space propulsion goals must be produced by a tool that 
accounts for performance and multiple, spatially resolved environmental variables. 
Efficient, validated CFD codes that model sufficient injector physics are necessary to 
meet these requirements. 

Generation of large amounts of complex information by these codes produces the 
need for a means to manage the data. The injector designer must be able to confidently 
and efficiently sort through this database to locate an acceptable compromise design. 
Global sensitivity and approximation tools can guide the designer effectively. 

Grid Sensitivity Investigation 

Initially the 54 cases identified by DOE were computed on an axi-symmetric 
geometry with 336x81 nodes. Only 33 out of the 40 fitting cases gave valid results. 
Results of the remaining seven cases contained unsteady features, which do not represent 
solutions of the steady state model employed. The RSAs generated with these data for 
TT max and X cc had R a 2 values of 0.961 and 0.976, respectively, suggesting a less than 
desirable fit. On checking the grid distribution in the combustion zone for the 33 cases 
used, it was determined that the grid resolution was insufficient. After a series of tests 
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involving addition of grid points in the axial direction in the combustion zone, a 430x81 
grid was found to be appropriate and used for the second run of the optimization study. 

To highlight the grid refinement, a comparison of the grid distributions is shown in 
Figure 4.8. The final grid was the product of tripling the axial node density in the 
combustion zone. The thick lines show the initial grid density, while the thin lines show 
final grid density. Note that, for clarity, only every sixth j-line is shown. New RSAs were 
generated from new solutions obtained on the fine grid. The new fits for TTma X and Xc C 
had R a 2 values of 0.989 and 0.995, respectively, representing a considerable improvement 
over the RS A performance based on the first, coarser grid. This time, only two out of the 
40 designs failed to produce valid results. It is to be noted that the R a 2 values for TF max 
and TW 4 are 0.999 for both the initial and final grid. This experience indicates that in 
addition to facilitating design optimization, the RSM can also help address the adequacy 
of the CFD solution accuracy. It offers insight into potential problems, based on the 
statistical regressions, from which the computations can be refined, thus improving the 
fidelity of the individual and collective databases. While this approach does not guarantee 
universally satisfactory outcomes, it does suggest clear directions to assess the critical 
area of data quality. 



Figure 4.8: Comparing the unrefined (336x81) {thick lines} and refined (430x81) {thin 
lines} grids. 
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Response Surface Generation 

In the following, data for each design variable from the 38 acceptable fitting cases 
were used to generate the RSAs. Both full and reduced quadratic polynomials are 
generated. Table 4.5 identifies o a and cr for the four objective functions and different 
polynomials generated. The full quadratic model is consistent in performance in terms of 
both <7 a and a. The reduced quadratic models either have a poorer value of <xor off er on jy 
marginal improvement over the full response surface model. Since there is no appreciable 
improvement by reducing the fits, there may be noise in the data or the quadratic models 
(both full and reduced) do not sufficiently represent the data. 


Table 4.5: Performance of full and reduced quadratic RSAs for the four objective 
functions (non-normalized values used). 




Full quadratic 

Reduced quadratic 

TF max 

Number of 
observations 

38 

38 ^ 


CT a 

0.00566 

0.00546 


a (14 points) 

0.00460 

0.00463 


Mean 

0.495 

0.495 

tw 4 

Number of 
observations 

38 

38 


<*a 

0.00803 

0.00795 


ct (14 points) 

0.00669 

0.00799 


Mean 

0.514 

0.514 

TT max 

Number of 
observations 

38 

38 


a a 

0.0413 

0.0401 


a (14 points) 

0.0396 

0.0382 


Mean 

0.560 

0.560 

Xce 

Number of 
observations 

38 

38 


CTa 

0.0205 

0.0197 


o (14 points) 

0.0178 

0.0186 


Mean 

0.497 

0.497 


Comparing the full quadratic model predictions for the fitting cases to the CFD 
results of the various objectives, the variations for TF max , TW 4 and Xc C were found to be 
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negligible (Figures 4.9A, 4.9B and 4.9D), suggesting no need for further improvement. 
However, a large number of points lie away from the best fit in the plot of TT max (Figure 
4.9C). 




Figure 4.9: Comparison between the best fit possible and as predicted by quadratic 

response surface. Optimum refers to the case when RSA and CFD values are 
identical. RS-CFD represents the value as for the current case (normalized 
Values Shown). A) TFmax- B) TW 4 . C) TTmax- D) Xcc- 

Based on this observation, a cubic model is generated for TT max . A full cubic fit in 
a 4 -design variable model requires a minimum of 35 designs points. To obtain a good fit, 
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the number of design points should be considerably larger than the required number. 
Since there are only 38 design points available from the fitting set, the testing set is also 
included in the fitting set. The PRESSrms is used to estimate the performance of the 
generated RSA along with the o a . Using the 52 design points, a full cubic was generated. 
The values of a a and PRESSrms were 0.0348 and 0.0598, respectively. The difference 
between the two measures of error is noticeable. Hence, a reduced cubic model is 
generated with the available design points. This improves the fit with values for cr a and 
PRESSrms being 0.0303 and 0.0388, respectively. Table 4.6 compares the performance of 
the full quadratic and the reduced cubic models for TT max . The reduced cubic model is 
seen to perform better statistically than the quadratic model. Hence, this cubic model is 
used for TT ma x in the optimization studies that follow. The RSAs for all four objectives 
are shown as Eqs. (4. l)-(4.4). These RSAs can then be used for sensitivity and 
optimization study. The RSAs are based on the scaled values of the design variables and 
objectives. 

Table 4.6: Performance of RSAs for the TT ma x- Reduced cubic RSA has 21 coefficients 


(non-normalized values used). 




Full quadratic 

Reduced cubic 

TT max 

Number of observations 

38 

52 



0.0413 

0.0303 


PRESSES 

0.0521 

0.0388 


Mean 

0.560 

0.591 


TFmax = 0.692 + 0.477(a) - 0.687(AHA) - 0.080(AOA) - 0.0650(OPTT) - 0. 167(a) 2 - 
0.0129(AHA)(a) + 0.0796(AHA) 2 - 0.0634(AOA)(a) - 0.0257(AOA)(AHA) + 


0.0877(AOA) 2 - 0.052 l(OPTT)(a) + 0.00156(OPTT)(AHA) + 0.00198(OPTT)(AOA) + 
0.0184(OPTT) 2 . (4.1) 
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TW 4 = 0.758 + 0.358(a) - 0.807(AHA) + 0.0925(AOA) - 0.0468(OPTT) - 0.172(a) 2 + 
0.0106(AHA)(a) + 0.0697(AHA) 2 - 0.146(AOA)(a) - 0.0416(AOA)(AHA) + 
0.102(AOA) 2 - 0.0694(OPTT)(a) - 0.00503(OPTT)(AHA) + 0.015 l(OPTT)(AOA) + 
0.0173(OPTT) 2 . (4.2) 

TT max = 0.370 - 0.205(a) + 0.0307(AHA) + 0.108(AOA) + 1.019(OPTT) - 0.135(a) 2 + 
0.0141(AHA)(a) + 0.0998(AHA) 2 + 0.208(AOA)(a) - 0.0301(AOA)(AHA) - 
0.226(AOA) 2 + 0.353(OPTT)(a) - 0.0497(OPTT)(AOA) - 0.423(OPTT) 2 + 
0.202(AHA)(a) 2 - 0.281 (AO A)(a) 2 - 0.342(AHA) 2 (a) - 0.245(AHA) 2 (AOA) + 

0.281 (AO A) 2 (AH A) - 0.184(OPTT) 2 (a) - 0.281(AHA)( a)(AOA) (4.3) 

Xcc = 0.153 - 0.322(a) + 0.396(AHA) + 0.424(AOA) + 0.0226(OPTT) + 0.175(a) 2 + 
0.0185(AHA)(a) - 0.0701(AHA) 2 - 0.251(AOA)(a) + 0.179(AOA)(AHA) + 
0.0150(AOA) 2 + 0.0134(OPTT)(a) + 0.0296(OPTT)(AHA) + 0.0752(OPTT)(AOA) + 
0.0192(OPTT) 2 (4.4) 

Optimization Process 

The generated RSA are used to conduct an optimization exercise and to study the 
relationship between the design variables and the objectives that are indicators of 
life/survivability and performance. Also, the ability to accommodate design features that 
promote extended life/survivability in a design while maintaining reasonable 
performance is explored. Three separate optimization studies, using the approach of 
desirability functions as explained in chapter 3, are presented below. First, single- 
objective minimizations are shown. Secondly, multi-objective optimizations are 
performed with equal weights. Finally, the multi-objective optimizations are conducted 
with variable weights. The obtained optimum solutions are compared with CFD 
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computations. Since the objectives are scaled to 0(1), the error measures have to be 
scaled accordingly to estimate the accuracy of the obtained solutions. In terms of the 
actual values the error for an objective is defined as 


error = 


_ ! yeFD }’RS 1 
y CFD 


(4.5) 


where y C FD is the solution obtained from the CFD analyses and )>rs is the prediction of the 
RSA. Using simple mathematics, not shown here, the error in the scaled variables can be 
written as 


error = 


y cfd y rs 

y CFD K 


(4.6) 


where the bar represents the scaled values, and K is defined as 


K = — ^ — (4.7) 

y max y min 

Here y m i„ and y max are the actual minimum and maximum values, respectively, 
based on the available set of fitting and testing data for that objective from the CFD 
analyses. In the dissertation, no bar is used for the notations used to represent the scaled 
design variables and objectives to avoid any confusion. Non-normalized values of the 
objectives refer to the scaled values and the normalized values refer to the re-scaled 
values of the objectives based on the maximum and minimum of the RSAs. This is to 
make some of the plots in the result section more meaningful. 

Single-objective optimization 

The purpose of this effort is twofold. First, the goal is to verify the performance of 
the optimization methodology in locating the minimum values for single objectives in the 
chosen design space. This verification is straightforward and necessary, but not sufficient 
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to conclude that the technique is useful for injector design. Table 4.7 shows the results of 
optimizing (here, minimizing) each objective separately. The numbers in parentheses in 
Table 4.7 indicate the weights applied during the optimization process. Here, a weight of 
one means that objective was included while a weight of zero means that objective was 
excluded. Based on the normalized values of the objectives the minimum value for each 


is necessarily 0. The optimizer found the correct value for each of the four cases. 


Table 4.7: Minimizing individual objectives. Value in parenthesis (1) indicates which 



Secondly, this process is helpful in understanding the injector operational 


influences via trend identification and variable groupings. The results from Table 4.7 are 


shown graphically in Figure 4. 10 to facilitate the discussion. It was shown earlier that the 
flow entrained in the large recirculation zone (Figure 4.7) regulated both the maximum 
temperature on the injector face, TF max and the chamber wall temperature, TW4. The 
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results from Opt-Cases 1 and 2 support this conclusion. Minimization of TF max very 
nearly minimizes TW 4 and vice versa. Not surprisingly, the designs for the two cases are 
also quite similar. Both designs are shear coaxial (a = 0) with the hydrogen flow area, 
AHA, and the oxidizer post tip thickness, OPTT, at their maximum values. Table 4.7 
shows that AOA is the only inconsistent variable in the two designs. In terms of the other 
objectives, when either TF max or TW4 is minimized, the maximum oxidizer post tip 
temperature, TT max , is high. The resulting moderate-to-long combustion lengths, Xc C , are 
consistent with the relatively slow mixing expected from a shear coaxial element. 
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Figure 4.10: Minimization of different objectives. (Case 1: TF max , Case 2: TW4, Case 3: 
TT max , Case 4: Xc C . Normalized values shown). A) Objectives. B) Design 
variables. 

Reference to Opt-Cases 3 and 4 in Table 4.7 and Figure 4.10 also suggests a 
correlation between TT max and Xc C , although not as tight as the one for TF^x and TW4. 
When TT max is minimized, Xc C is low and vice versa. Further investigation needs to be 
done to completely understand the physics that underlies this correlation. Here, both 
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designs are impinging-like with the hydrogen flow angle, a, at or near its maximum 
value and AHA and OPTT at their minimum values in the chosen design space. Again, 
Table 4.7 shows that AOA is the only inconsistent variable between the two cases. Both 
minimizations result in very high values for TF max and TW 4 . Not surprisingly, the 
impinging elements represented by the two designs yield significantly shorter combustion 
lengths than the shear coaxial designs. 

The trends for a, AHA and OPTT are consistent between the two pairs of 
dependent variables but AOA varies among the four cases. With the other design 
variables set at the optimum design levels, Figure 4.1 1 A shows the variation of TF max and 
TW 4 with AOA. A similar plot for TT max and Xc C is shown in Figure 4.1 IB. The 
maximum injector face temperature, TF max , is least sensitive to AOA. It is also the only 
objective to exhibit a minimum value in the interior of the design space. This finding 
suggests that expanding the design space could result in more robust designs for the 
multi-objective optimization. The trends for TW 4 and Xc C are similar to each other, with 
AOA zero for both cases. The trend for TT max relative to AOA is opposite. 

Table 4.7 and Figure 4. 10 show the resulting design is different for each of the four 
single-objective optimizations. These results indicate that designing robust injectors 
demands consideration of all the independent variables together. Figure 4.10 clearly 
shows the competing design trends. Thus, meeting a set of design requirements will 
require a compromise design. 
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Figure 4.1 1: Variation of objectives with respect to AOA for fixed a, AHA and OPtt 

(normalized values shown and D is the desirability function. All objectives are 
equally weighed in the composite function). A) TF max and TW 4 for a= 0 , 
AHA=1 and OPTT=l. B) TT max and X« for <x=l, AHA=0 and OPTT=0. 

Multi-objective optimization 

To concurrently evaluate component life/survivability and performance 
considerations, a multi-objective optimization study is carried out. Starting with 
performance, X cc , as most of the injector design tools do, the objectives influencing 
thermal environments are added one at a time to study the effect on the resulting 
optimum designs. These optimum designs are presented in Table 4.8 and Figure 4.12 
where Opt-Case 4, with Xc C minimized is repeated as the starting point. The design for 
this case is an impinging element (a = 0.917) with minimum flow areas and the thinnest 
oxidizer post tip. The consequences of this design are a minimum Xc C , a low TT max and 
very high values of TF max and TW4. 

When TT m ax is minimized along with X cc in Opt-Case 5, the optimum design (foes 
not change significantly from that in Opt-Case 4. The hydrogen flow angle increases 
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form 0.917 to 1.0. The values of AHA, AOA and OPTT remain unchanged from Opt- 
Case 4. The combustion length, X cc , is unchanged as compared to Opt-Case 4, whereas 
TT max is marginally improved from 0.181 to 0.152. The other thermal objectives remain 
at very high levels. Recall from Table 4.7 that AOA was 1.0 for the single-objective 


optimization of TTmax. 

Table 4.8: Study of effect of life/survivability on performance. Value in parenthesis 



a 

AHA 

AOA 





Xcc 

4 

0.917 

0 



mu 

mtm 

0.182 

(0) 


CFD 




• 

0.981 

eeqh 


-0.004 

iSH 





0.08 



0.12 

5 

1 

0 

0 

0 

1.0 

(0) 


mm 

mi 

CFD 





0.998 

0.927 


-0.004 






0.03 

0.01 


0.14 

6 

1 

1 

0 

0 

0.376 

(1) 

0.224 

JO) 

0.155 
0) 

mm 

CFD 





0.369 

0.214 

0.252 

0.264 









0.38 

i 

1 

1 

0 




IU 

mm 

CFD 





0.369 

0.214 

0.252 

0.264 

Error 

(%) 





0.10 

0.12 

3.93 

0.38 


More insight regarding the optimum design can be gained by visualizing Figure 
4.1 IB. Although the individual TT max and X cc drive the design towards the opposite ends 
of the range of AOA, the composite desirability function, D, has marginal variation with 
respect to AOA. The optimizer picks the value of AOA that gives the maximum value of 
the desirability function, D. The optimum design at this point is affected mostly by the 
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other design variables. This shows the benefits of using a global optimization procedure, 
which helps identify the trend of the composite function over the design space even when 
the individual objectives have opposing trends. 
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Figure 4.12: Composite minimization of objectives with different weightings. (Case 4: 
(0,0,0, 1), Case 5: (0,0, 1,1), Case 6: (1,0, 1,1), Case 7: (1,1, 1,1). The values in 
parenthesis indicate weights for (TF max , TW 4 , TT max , X cc ). Normalized values 
shown). A) Objectives. B) Design variables. 

Opt-Case 6 minimizes TF max simultaneously with TT max and X cc . In terms of the 
design variables, AHA has now shifted from 0 to 1, while the other three remain at the 
same values as in Opt-Case 5. This shift in AHA is consistent with the single-objective 
minimization of TF max shown as Opt-Case 1 in Table 4.7. This single change in the 
design produces dramatic changes in the objectives. As desired, the value of TF max has 
decreased significantly from 1.0 to 0.376. The concurrent large drop in the value of TW 4 
is consistent with the earlier conclusion linking TFmax and TW 4 . The oxygen post tip 
temperature, TT max , remains essentially unchanged. The increase in X cc from 0 to 0.279 is 
a negative impact of the design change. Figure 4.13 sheds more light on the situation. 
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Figure 4. 13 A plots the individual objectives as a function of AHA. Figure 4.13B is 
similar, with the individual objectives plotted as a function of AOA. These figures clearly 
illustrate that both TF max and TW 4 are considerably more sensitive to AHA than AOA. 
Accordingly, the optimizer uses the hydrogen flow area to most efficiently regulate 
TF max . Figure 4.13A also shows the positive slope of Xc C relative to AHA, which is 
responsible for the fairly large performance drop seen in this design. 




Figure 4. 13: Variation of objectives with respect to AHA or AOA with other design 

variables fixed, (normalized values shown and D is the desirability function. 
All objectives are equally weighed in the composite function). A) TF max , TW4, 
TT max and Xc C , with respect to AHA for a=l, AOA=0 and OPTT=0. B) TF max , 
TW4, TT max and Xc C , with respect to AOA for a=l, AHA=1 and OPTT=0 

Opt-Case 7 adds the final objective, TW4. Table 4.8 shows that addition of TW4 has 
no effect on the design and consequently no effect on any of the objectives. The chamber 
wall temperature, TW 4 , has previously been shown to be linked with TF raax , so this result 
is not surprising. 

The injector design resulting from inclusion of all the objectives in the optimization 
process (Opt-Case 7) is different from Opt-Case 4 where only the performance indicator, 
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Xcc was optimized. In the design space evaluated for this effort, AHA has a very strong 
effect on the system, especially for TF max and TW4. Reference to Figure 4.13 shows that 
at the zero value for both AHA and AO A, X cc is still decreasing. Performance could 
probably be further improved (i.e., Xc C shortened) by enlarging the design space to 
include smaller values of AHA and AOA. In the process of designing a real injector, the 
large drop in performance affected by the consideration of TF max and TW4 would 
probably require a compromise design with certain variables weighted more heavily than 
others. 

To further probe the interplay between injector performance and component 
life/survivability, a composite optimization study is conducted with varying weights on 
the individual objective functions embodied by the composite desirability function. This 
is a functionality that the new design optimization methodology must possess. The 
designer must be able to weight, or favor, one or more dependent variables to have 
maximum flexibility in addressing design requirements. The use of geometric mean and 
desirability functions is an approach which addresses such a requirement. Additionally 
MOGA have also been shown to work well for such problems. The results of the MOGA 
will be presented in the later sections of this chapter. The results of the current study are 
shown in Table 4.9 and Figure 4.14. Reference to Table 4.9 (the numbers in parenthesis 
indicate the weight in the composite desirability function for that variable) shows that the 
designs suggested by Opt-Case 8 and Opt-Case 10 are similar. In Opt-Case 8 the 
life/survivability influencing objectives TF max , TW 4 , and TT max are weighted by a 2:1 
ratio over the performance indicator (X^). Opt-Case 10 is the opposite, with performance 
having a 2:1 weighting over the life/survivability objectives. Reference to Tables 4.8 and 



91 


4.9 shows that the design variables and objectives for these two cases are almost identical 
to those in Opt-Case 7, where all four objectives are weighted equally. Figure 4. 13B 
shows the composite desirability function D, for Opt-Case 7, to be reasonably flat as a 
function of AOA. Thus, small weightings have little or no effect on the injector design. 
Table 4.9: Study of influence of life/survivability and performance objectives on each 


other. Value in parenthesis indicates the weighting on the objective functions 
(normalized values shown). 


Opt- 

Case 

a 

AHA 

AOA 




TTmax 

x cc 

8 

i 

1 

0.192 

0 

0.346 

(1) 

0.210 

(1) 

m 

E 

niEfli 

CFD 





0.343 



0.326 

Error 

(%) 



|B 


0.04 



0.22 

9 

0 

1 

0.497 


II 

PH 

0.471 

.(5) 

E 

EH 

CFD 





msm 


0.404 

0.613 

..Bnt 






0.07 

2.53 

0.32 

10 

1 

1 

0 

0 


ICB «■ 

IJEBI 


CFD 





0.369 


0.252 

0.264 

13W 

1531 1 





0.10 

0.12 

3.93 

0.38 

ii 

0.612 

0 

0 




0.282 

(0.1) 

EH 

CFD 





0.911 

0.890 

0.361 

0.0070 

!■ 





0.10 

0.09 

3.10 

0.16 


The other two cases shown in Table 4.9, Opt-Case 9 and Opt-Case 11, have relative 
weightings of 50:1 and 1:50, respectively for life/survivability (TF raax , TW 4 , and TT max ) 
and performance indicators (X cc ). For Opt-Case 9, the design tends toward the cases 
where TF ma x and TW 4 were individually minimized. The exception is that OPTT is equal 
to zero, which is required to minimize TT max . For Opt-Case 1 1, the resulting design tends 
toward Opt-Case 4 where Xc C is minimized. Even with this strong weighting, the effect of 
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including TF raa x and TW4 is felt with the value of a equal to 0.612. When X ec is 


minimized by itself, a is equal to 0.917. 
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Figure 4.14: Composite minimization of objectives with different weightings. (Case 8: 
(1,1, 1,0.5), Case 9: (5, 5, 5, 0.1), Case 10: (0.5, 0.5, 0.5,1), Case 11: 

(0.1,0. 1,0. 1,5). The values in parenthesis indicate weights for (TF max , TW4, 
TT max , Xc C . normalized values shown). A) Objectives. B) Design variables. 

Figure 4.15 shows the joint desirability function, D, for Opt-Case 9 (with a 
weighting of 50: 1 in favor of life/survivability over performance), to be fairly flat. 
Accordingly, small weightings have little effect on the design. Significant weightings 
must be applied to push the design very far toward either life/survivability or 
performance. Opt-Case 1 1 has good performance and a low injector tip temperature. 
However, the chamber wall and injector face temperatures are very high. In an actual 
design, if this level of performance was required, active cooling could be required for the 
chamber wall and injector face. 

Additional CFD solutions were executed to confirm the optimum designs obtained 
from RSA for all 1 1 Opt-Cases. These results are shown in Tables 4.7, 4.8 and 4.9. The 
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objective error, as defined by Eq. (4.5), is also shown for each case. The error is on the 
order of four percent or less in the case of TT mBX . The error for the other three objectives 
is less than one percent. The larger discrepancy in TT max is likely due to the steady state 
assumption made in the CFD analyses for this effort. The injector problem is actually 
unsteady. 


Response Functions 



Figure 4.15: Variation of TF max , TW 4 , TT max and Xc C , with respect to AOA for a=0, 
AHA=1 and OPTT=0 (normalized values shown and D is the desirability 
function. Case 9: Temperatures weighted over performance in the ratio of 
50:1). 

Preliminary Observations 

Based on the preliminary optimization studies specific observations can be made. 
Single-objective optimization. The optimizer reliably locates single objective 
minimum values. Minimization of the individual objectives yields a different injector 
design. This establishes the fact that designing an efficient long life/survivability injector 
requires concurrent evaluation of all the relevant design variables. 

For the current study, 
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1 . The four objectives, based on the design obtained during their individual 
minimization, largely fall into two groups: (TF max , TW 4 ) and (TT max , Xc C ). 
However, there are differences between them. In each group, the responses of 
the two design objectives to AOA are different. The observation also 
indicates that the three life/survivability-related objectives require 
compromises between design variables. 

2. Minimizing TF max and TW 4 leads to a design with a equal to zero (shear 
coaxial element), maximum fuel flow area and thickest post tip. This design 
also yields moderate to poor performance due to the slow mixing across the 
shear layer. 

3. Minimizing TT max and X cc results in an impinging-like design with a equal to 
one. It also has the minimum fuel flow area and the thinnest post tip 
thickness. This design performs well, but has very high wall and injector face 
temperatures. 

Multi-objective optimization. The injector design, when multiple objectives are 
included, is different from any of the single objective minimizations. For example, 
although the individual objectives, TT max and X cc drive the design towards the opposite 
ends of the range of AOA, the composite desirability function accounting for all their 
effects exhibits only marginal variation with respect to it. The optimum design is 
affected mostly by other design variables. 

One can clearly see the benefit of using a global optimization procedure, which 
helps identify and interpret the trend of the composite function over the design space, 
especially when the individual objectives have opposing patterns. 

For the current study the following specific observations can be made. 

• Equal weights 

1 . The design with all 4 objectives included (Opt-Case 7) is still an impinging 
element with low to moderate temperatures, but marginal performance. 

2. Figure 4. 13 indicates that lowering either propellant flow area beyond the 
current limit would probably decrease Xc C and thus increase performance. 
Lowering the oxidizer area would probably have less adverse effect on TF max 
and TW 4 . 
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• Unequal weights 

1 . With modest weights on either the performance (Xc c ) or life/survivability 
(TFnmx, TW 4 and TT max ) the design does not change appreciably. 

2. Opt-Case 1 1 (with a large emphasis on performance) gives very good 
performance, modest TT max , but very high values of TT max and TW 4 . If 
enlarging the design space does not help, active cooling may be required. 

Sensitivity and Trade-Off Analysis 

The analysis is divided into two parts. Firstly, the global analysis is carried out 
where the entire design space is explored to estimate the sensitivity of the objectives to 
the design variables. Following this the POF is explored to understand the trade-offs 
between the objectives and also the different design trends. 

Global analysis 

The global analyses addresses, both the sensitivity of objective variability to design 
variables and the interaction between objectives over the complete design space. The 
global sensitivity indices are computed using the variance-based non-parametric 
approach proposed by Sobol (1993) and described in chapter 3. The code written using 
Matlab (2002) was validated using a well-known benchmark problem (Mckay, 1997). 
Table 4. 10 summarizes the results of this study and lists the essential and non-essential 
design variables with respect to individual objective variability. A design variable is 
considered essential if it is responsible for at least 5% of the objective variability. 

Figure 4.16 shows the percentage of main factor (Si) contribution of different 
design variables to individual objective variabilities. The variability of TF max is largely 
influenced by AHA and moderately by a (Figure 4.16A). The effect of the other design 
variables is marginal suggesting that they are non-essential and in principle could be 
fixed. The variability of TW 4 is considerably influenced only by AHA (Figure 4.16B). 
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TT max is influenced considerably by OPTT and marginally by a (Figure 4. 16C). For X cc , 
AHA, AOA and a have considerable influence (Figure 4.16D). The total sensitivity 
indices ( S ‘ otal ) were computed and compared with the main factors (Si) and it was found 
that the contributions of the cross-interactions among the design variables to the 
objectives variability were negligible. 

Table 4.10: List of essential ( V) and non-essential (x) variables for each objective and the 


mean errors between the modified RSA and original RSA. An essential 
variable accounts for at least 5% of the objective’s variability. 


Objective 


Design 

a 

Variables 

AHA 

AOA 

OPTT 

Mean 
Error (%) 

Life/Survivability 

TF max 

V 

V 

X 

X 

5.2 

indicators 

tw 4 

X 

V 

X 

X 

11.5 


TTmax 

V 

X 

X 

V 

6.7 

Performance 

indicator 

Xc C 

V 

V 

V 

X 

6.1 


TFmax TW4 


OPTT, 



AHA, 0.955 


A B 

Figure 4.16: Main factor (Si,) influence on objective variability. A) TF max . B) TW 4 . C) 
TT m ax- D) X cc . 
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TTmax 


Xcc 



Figure 4.16: Continued. 

This global sensitivity analysis shows that different design variables can have 
varied effect on individual objectives. Such a study can help the designer fix some of the 
design variables while analyzing the design. For example, in the current injector design 
study, for objective TF max , the design variables AOA and OPTT can be fixed at their 
mean value (0.5) as this do not result in significant differences in the prediction. The 
mean error (difference between predictions) throughout the design space between 
modified RSAs (obtained by fixing the non-essential design variables at their mean 
values), and the RSAs listed in Eqs.(4.1)-(4.4), (referred to as original RSA in the rest of 
the text) are given in Table 4.10. The error for TW 4 is about 12% suggesting that the cut- 
off for the non-essential variables may have to be lowered to capture additional features 
of the original model. The mean error in the modified RSA of the remaining objectives is 
about 6%, suggesting that they capture the original RSA reasonably well. This 
information can be used for dimensionality reduction and therefore, to ease the search for 
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optimum designs. Similar concept has been addressed in the past using a different 
approach. For example, Knill et al. (1999) have used linear aerodynamic to identify the 
important terms in the polynomial-based RSA which were then used for creating the 
surrogate model from Euler analyses. This reduced considerably the number of points 
needed for the Euler design of experiments. For the current work since the number of 
design variables is small the remaining studies are carried out without fixing the non- 
essential variables. From an engineer’s viewpoint and interest, this study provides an 
insight into the physics of a design problem by highlighting the design features that 
govern the individual objectives. 

Linear surrogate models for the four objectives are constructed as a function of the 
four design variables. The corresponding coefficients are shown in Table 4. 1 1 . The 
magnitude of the coefficients agrees well with the essential and non-essential nature of 
the design variables for each objective. Additionally, the R 2 values for the linear RSA are 
compared with those of the original RSAs (referred as R 2 nonlinear in Table 4.1 1). 
Comparing R 2 unear with R 2 nonlinear shows that most of the variability is accounted for by 
the linear model and the additional terms in the nonlinear model give marginal 
improvement. 

Table 4.1 1 : Coefficient associated with the different terms in the linear RSA and 


comparison ofi? 2 values of linear RSA with that of nonlinear RSA. 



a 

AHA 

AOA 

OPTT 

Intercept 

R linear 

p2 

nonlinear 

TFmax 

0.237 

-0.634 

-0.033 

-0.066 

0.733 

0.990 

0.999 

tw 4 

0.0719 

-0.775 

0.114 

-0.052 

0.826 

0.986 

0.999 

TT max 

-0.222 

0.108 

-0.063 

0.662 

0.367 

0.918 

0.994 

Xcc 

-0.234 

-0.402 

0.433 

0.108 

0.138 

0.958 

0.997 
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A correlation analysis was then carried out by Goel et al. (2004) to observe the 
interactions between objectives. A total of 14641 (1 1 4 ) design points were generated over 
the complete design space by varying one variable at a time by a constant value. The 
objectives were calculated using the original RSAs. The correlation matrix Cdes and 
corresponding p- values are computed using Matlab. 

The correlation matrix, Cdes shows that there is a strong correlation between 
objectives TF max and TW4, as the corresponding coefficient is very close to 1. P-values 
and 95% confidence intervals for the correlation coefficients also establish the statistical 
significance of the results. Low P-values (« 0.05) confirms the significance of the 
correlation results. This finding is in agreement with the observations made in the 
preliminary optimization study. The combustion chamber wall temperature, TW4, is 
excluded and the optimization problem is formulated with the remaining three objectives. 

' TF max X cc TW 4 TT^ 

TF max 1.00 -0.773 0.947 -0.272 
Cdes = x„ -0.773 1.00 -0.583 0.263 
TW 4 0.947 -0.583 1.00 -0.193 
TT -0.272 0.263 -0.193 1.00 

max 

Pareto front analyses 

The trade-offs between objectives and sensitivity analyses are carried out at the 

POF. The Pareto fronts of (TF max , Xc C ) and (TF max , TTmax) are of interest as these 
objectives conflict one another. Goel et al. (2004) have generated the POFs with the aid 
of NSGA-II. Figure 4.17A shows the relation between TFmax and Xc C . It can be seen that 
the POF in this case is linear over a large region. A small increase in the value of TF max 
(« 10%) reduces the combustion length, X cc , by nearly 50%. Figure 4.17B shows the 
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relation between TF max and TT max . It is obvious that the POF is non-convex. It is also 
seen that a small drop in the value of the face temperature (« 10%) can reduce the tip 
temperature TT max by nearly 60%. Hence at a small cost of TF max both X cc and TT max 
improves considerably. 



A B 


Figure 4.17: Two-objective Pareto optimal front. A) TF max vs Xc C . B) TF max vs TT max . 

Following the trade-off studies, the three-objective (TFmax, Xc C and TT max ) 
optimization problem is solved using the NSGA-II algorithm (Goel et al., 2004). The 
solutions obtained are further refined using the e-constraint local search strategy coded 
in Matlab. In this strategy, first, one of the objectives TF max , is treated as the objective 
function and rest of the objectives X« c and TT max are treated as equality constraints. The 
constraint value is set equal to the corresponding objective value as found by NSGA-II 
simulation. The corresponding design variable vector is used as initial guess. This 
procedure is repeated for all individuals in the population. This gives a set of Pareto 
optimal solutions, referred to as Set A in this study. Next, Xc C is chosen as the objective 
and the other objective functions TT max and TF max are treated as constraints. For this 
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problem, the obtained Pareto optimal set is referred to as Set B. Similarly, the Pareto 
optimal set, Set C, is obtained with TT max as objective and TF max and X cc as constraints. 
Now all the Pareto optimal solutions obtained so far; i.e.. Sets A, B, C and the original 
NSGA-II set are combined. To find the true Pareto optimal front, a non-domination check 
is carried out on this set of 400 solutions. This yields 254 non-dominated solutions. After 
removing the duplicates, there are 249 solutions, which are on the POF. These solutions 
are shown in Figure 4.18. This solution set is the global Pareto optimal solution set. The 
optimal solutions listed in Tables 4.7-4.9 are also plotted in Figure 4.18 



Figure 4.18: Pareto optimal solution set (+) and the multi-objective optima listed in 
Tables 4.7-4.9 (o). 

The optimal solutions shown in Tables 4. 7-4.9 lie along the POF obtained using 
MOGA and the local search algorithm. Although identical solutions are not identified by 
the different multi-objective optimization approaches, the trend of the Pareto front is 
adequately captured. This shows the effectiveness of the different approaches. 
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Goel et al. (2004) has then used the hierarchical clustering algorithm (Jain and 
Dubes, 1988) to divide the obtained POF into 9 clusters for sensitivity and trade-off 
analyses. Values of the design variables and objectives for these designs are shown in 
Table 4.12. Graphically, the solutions are shown on the POF in Figure 4.19. It is clear 
from Figure 4.19 that solutions are selected uniformly over the design space. In Figure 
4.20, the box plots of the design variables for clusters 1, 3, 6 and 9 are shown. The box 
plots for the objectives are shown in Figure 4.21. These plots highlight the variability of 
the design variables and objectives in each cluster. 


Table 4.12: Objective function and design variables of nine (9) representative designs 


from the Pareto optimal solution set. 

Cluster 

a 

AHA 

AOA 

OPTT 

TFmax 

Xc C 

TTmax 

1 

0.000 

1.000 

0.842 

0.712 

0.023 

MEM 

0.880 

2 

0.000 

1.000 

0.356 

0.587 

0.028 

Mm 

0.0890 

3 

0.000 

1.000 

0.442 

0.014 

0.054 

wttm 

0.452 

4 

0.094 

1.000 

0.000 

0.015 

0.126 

0.453 

0.466 

5 

0.668 

1.000 

0.732 

0.000 

0.259 

0.681 

0.229 

6 

0.600 

0.670 

0.000 

0.000 

0.489 

0.264 

0.226 

7 

0.295 

0.108 

0.000 

0.354 

0.719 

0.129 

0.641 

8 

0.314 

0.066 

0.000 

0.055 

0.776 

0.097 

0.357 

9 

1.000 

0.014 

0.680 

0.000 

0.935 

0.138 

-0.043 


For cluster 1 it is seen that the value of a is fixed at 0 (shear coaxial injector) 
(Figure 4.20A) and AHA is fixed at 1 (Figure 4.20B). This suggests that in this cluster, 
the designs are sensitive to a and AHA both of which reach their extreme values. The 
remaining two design variables AOA (Figure 4.20C) and OPTT (Figure 4.20D), vary 
over a range. It is observed that TF ma x is minimized (Figure 4.21 A) where as Xc C and 
TT m ax lie near their maximum and have little variability. Hence the designs in cluster 1 
tend to minimize TF max and represent shear coaxial injector designs. The design variables 
AOA and OPTT do not influence TF max but affect the remaining objectives, Xc C and 
TT max . Partial correlation coefficients are estimated to obtain the relationship between the 
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variances of these design variables and objectives (Table 4.13). It is noticed that as AO A 
increases, Xcc increases (Rcoit = 1-00) and TT max decreases (Rcon = -0.638), thereby 
requiring a compromise in the design. As OPTT decreases both Xc C and TT max decrease 
(Rcoit = 1 -00 for both). 



Figure 4.19: Pareto optimal solution set and nine (9) representative solutions from the 
same set. The circles identify the representative of a specific cluster. 




Figure 4.20: Box plots for the design variables in clusters 1, 3, 6 and 9. A) a. B) AHA. C) 
AOA. D) OPTT. 
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Figure 4.20: Continued. 

Similar observations can be made for clusters 3, 6 and 9. Figures 4.21A-4.21C 
show that with increase in TF max , from cluster 1 to cluster 9, Xcc and TT max decrease 
(based on the median of the box plots). This highlights the trade-off between the 
objectives. Cluster 9 provides information about the opposing trend as to what was 
observed in cluster 1. As objectives Xc C and TT max are minimized (Figures 4.21B and 
4.2 1C) an impinging injector design is obtained (a ~ 1, Figure 4.20 A) with TF max 
exhibiting high values (Figure 4.21 A). The AHA is near minimum (Figure 4.20B) 
contrary to the design in cluster 1 where high AHA minimized TF max . The AOA has 
considerable variation which suggests that the variability of objectives, X cc and TT max are 
not largely affected by this design variable. Figure 4.20D shows that Xc C and TT max are 
minimized for the minimum value of OPTT. 

Table 4.13 gives the partial correlation coefficients for the set of design variables 
and the objectives in each cluster which shows considerable variation. The partial 
correlation coefficients for the combinations left out are effectively zero. These measures 
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can be used to tune the design variables in a chosen cluster so as to improve on the 


objectives as per design requirements. 





Figure 4.21: Box plots for the objectives in clusters 1, 3, 6 and 9. A) TF max . 


TT 


max- 


B)Xc C . C) 


In this study an elaborate optimization study with different compromised design 


goals has been presented along with both a global and POF sensitivity analyses. The 
conclusion drawn from these studies will be presented in chapter 5. 
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Table 4.13: Partial correlation coefficients (Rcorr) of design variables vs. objectives for 


different clusters along the POF 

Cluster 

Design variable 

TF m ax 

x cc 

TT m ax 

1 

AOA 

- 

1.000 

-0.638 

1 

OPTT 

- 

1.000 

0.991 

3 

AOA 

- 

1.000 

- 

6 

a 

0.982 

-0.735 

-0.983 

6 

AHA 

-0.999 

0.994 

- 0.729 

9 

a 

0.877 

-0.203 

- 0.769 

9 

AHA 

-0.992 

0.983 

0.816 

9 

AOA 

-0.977 

0.997 

-0.940 


CHAPTER 5 

SUMMARY, CONCLUSION AND FUTURE WORK 

In this chapter the work done is this dissertation is summarized and conclusions 
drawn. Firstly, the summary and conclusion of the work related to Navier-Stokes (NS) 
CFD code verification is presented. Following this the design optimization study is 
summarized and conclusions drawn. Finally possible future research topics are briefly 
mentioned. 

Navier-Stokes Code Verification 

In this study, least square extrapolation (LSE) method has been mainly 
implemented on Navier-Stokes computations solved in the finite volume (FV) 
formulation. Two coarse grid solutions are extrapolated onto a fine grid using constant 
weighting parameters estimated by minimizing the residuals of the NS equations on the 
fine grid in the least square sense. The lid-driven cavity flow with two different boundary 
conditions (constant and variable lid velocities) are used to study the robustness and 
efficiency of LSE in situations of singularities, non-linearity and coupled equations. The 
initial test for LSE is done using a 2-D turning point problem is tested. Following this, 
only pressure is extrapolated for the lid-driven cavity flows. Finally LSE is implemented 
on the complete NS equations to extrapolate the velocity components and pressure 
simultaneously. 

The key observations pertaining to the LSE approach can be summarized as 
follows: 
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• It has been noticed with the aid of a 2-D turning point problem and a laminar lid- 
driven cavity flow that this method can yield solution improvement for linear as 
well as nonlinear PDEs. Additionally it has been shown with this method that the 
extrapolating weights based on minimization of the residual of the governing 
equation performs better than the weights used by RE based on assumed order of 
convergence. 

• It has been shown with the aid of two different Reynolds number (Re = 5.33 and 
1000) that LSE works well over a wide range provided the base grids captures the 
flow features adequately. 

• Minimizing the residuals does not always lead to minimizing the extrapolated 
solution error. In particular, it has been shown that flow with singularities Would 
introduce large variation in the eigenvalues of the set of residual equation that is 
used for least square approach. This would influence the accuracy of the solution 
even though least square might effectively reduce the L 2 norm of the residual. 

• In practice, neglecting the regions of singularity can improve the performance. This 
suggests the requirement of further study to address issues of singularities in 
different flows. 

• The LSE is seen to work well for the pressure-velocity coupled system. The Picard- 
like iteration adapted for the nonlinear momentum equation converges well f or the 
shown case. 

• It is also seen that the individual extrapolation of pressure or the coupled pressure- 
velocity extrapolation does not affect the accuracy of the extrapolated pressure field 
drastically. 

Overall issues related to singularities in the lid-driven cavity flow were noticed. 
This resulted in an elaborate study to understand the problems related to it and its 
influence on LSE. In this study only a constant weighting parameter is used. The 
preliminary implementation of spatially dependent weighting function exhibits 
inconsistent performance, which requires further investigation. 

Design Optimization Study 

In this study, the integration of CFD and optimization tools to explore the design of 
a single element rocket injector is presented. The design variables indicating the thermal 
environment and performance of the element are identified and the design points selected 
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using a DOE technique. CFD solutions are obtained for each design and appropriate 
orders of polynomial-based RSAs are generated for the individual objectives. Based on 
the performance of the RSAs a grid refinement study is carried out and a new, more 
efficient grid is generated. Utilizing the results obtained from CFD computations on the 
fine grid, a preliminary multi-objective optimization study is carried out using a 
composite function based on geometric mean approach and different design trends, 
correlations and trade-offs observed. 

Following the preliminary study, RSA, global sensitivity indices and genetic 
algorithms have been used to conduct an elaborate sensitivity and trade-off analyses of 
the injector design. Global analyses have been conducted to observe the influence of 
design variables on the objective variability and estimate the correlation between 
objectives. A Pareto optimal front (POF), generated by Goel et al. (2004) using a RSA- 
multi-objective genetic algorithm (MOGA) coupled approach, is used to conduct trade- 
off studies between objectives. A three-objective POF has been divided into 9 clusters 
and box plots have been used to observe the variability of the designs in each cluster. 
Additionally, partial correlation coefficients have been calculated to derive a relationship 
between the objectives and design variables in each cluster. Such analyses provide the 
designer with enough information to filter out designs based on his specific needs. Some 
observations based on the tools used, their performance, and the different design aspects 
of the injector are highlighted below. 

• The degree of spread of a dependent variable over a parametric space is an 
important criterion for validation experiment design. Evaluation of large parametric 
spaces provides insight into potential validation measurements. 

• Correlations of the variables highlight the trends and provide insight into physical 
processes that regulate the objectives. 
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• The study confirms that the details of the injector design govern the performance. 
Directing the fuel flow towards the element axis and thinning the oxidizer post lead 
to better performance. The same details are shown to have a large impact on the 
injector and combustion chamber environmental variables as well. 

• Using CFD to evaluate and include more design variables and objectives in the 
conceptual design phase has potential to create more robust designs. This increase 
in the dimensionality of the design problem also creates the requirement for an 
optimization tool to guide the designer through often opposing trends to an 
acceptable design compromise. 

• Based on the preliminary optimization study it was found that the proposed CFD- 
RSM approach worked well. It was fairly efficient in terms of the number of CFD 
solutions required and provided sensible results in view of the physics of the 
problem. 

• For both single- and multi-objective optimization efforts, in many cases, the 
resultant designs are at an extreme of the design space. This indicates the design 
space could be enlarged subject to practical design considerations. Again, a merit 
of the RSM is that the solutions obtained can be used repeatedly when revising the 
design scope. Furthermore, different scenarios, with either single- or multi- 
objective optimization, can make use of the information supplied by the response 
surface repeatedly. 

• The global sensitivity indices (main factor and total sensitivity indices) help 
identify which design variable is essential for objective variability. Additionally, 
point out the effect of cross interactions among the design variables on the 
objectives variability. 

• The variability of TF max depends on a and AHA, TW4 on AHA, TT max on a and 
OPTT and Xc C on a, AHA and AOA. The rest of the design variables for each 
objective can, in principle, be fixed, as their contribution to the objective variability 
is marginal. It is noticed that the mean errors between the modified RSAs and the 
original RSAs are about 6-12%. 

• It is observed that TF max and TW 4 are strongly correlated and hence TW 4 is 
excluded from the multi-objective optimization study. 

• The conflicting nature of the objectives (TF max , TT max ) and (TF max , Xc C ) is observed 
by examining the corresponding POFs. For a small cost in TF max (~10%) , both Xc C 
and TT max (~50%) can be decreased considerably. 

• The POF obtained for the 3-objective study (TFmax, TT max and Xc C ) gives an idea of 
the various compromises among the different objectives. A hierarchical clustering 
approach helps divide the POF into regions of similar objective goals. 
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• Box plots and partial correlation coefficients can assist the design selection along 
the POF. Box plots identify the variability of design variables and objectives in 
each cluster. Partial correlation coefficient identifies linear trends between the 
design variables and objectives variabilities. 

• To minimize TF ma x it is essential to have a shear coaxial injector (a = 0) and 
maximum change in baseline H 2 flow area. To minimize Xc C and TTmax, an 
impinging injector (a > 0) is required with minimum change in baseline H 2 flow 
area. 

T his study has offered in depth information about the different design aspects 
pertaining to the design variables and objectives. The importance of different design 
variables on individual objectives has been identified and the interaction between the 
objectives noticed. The POF of the 3-objective study gives different choices for the 
designer to look into but the key observation is that for a marginal cost in the injector 
face temperature, TF max , the performance, X cc and O 2 post tip temperature, TT max can be 
considerably reduced. 

The MDO tools, namely DOE and RSM, have been shown to work well along with 
CFD. Vaidyanathan et al. (2003b) have shown that these tools can be use to address 
model fidelity issues in the context of a turbulent cavitating flows. In this work they have 
used these tools to identify the optimum values of the empirical parameters present in the 
turbulence and the cavitation models. This highlights the broad scope of the tools 
presented in this dissertation. 

Future work 

Future research directions are briefly presented in this section. First issues related 
to LSE are addressed and possible strategies suggested. Next, ideas for effective use of 
design tools used in this dissertation are presented. 
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Navier-Stokes Code Verification 

In this dissertation, the weight function in LSE has not been allowed to vary 
spatially. For the problems tested and the techniques implemented, spatially dependent 
weight function was found to deteriorate the extrapolation and hence only a constant 
weight was used for all the study presented. One of the key issues during the 
extrapolation is the interpolation of coarse grid solutions onto the fine grid. This 
interpolation introduces high spurious wave number terms which combined with the 
spatially dependent weight function deteriorates LSE. For flows like lid-driven cavity 
flows, such high wave number components are amplified during residual estimation. The 
LSE approach is implemented on these polluted residuals and hence does not guarantee 
the minimization of the solution error. More details regarding this issue can be found in 
the work of Garbey and Shyy (2004). They have proposed using a low mode 
approximation of the residual for LSE implementation. This smoothing of the residual is 
shown to be effective in improving the performance of LSE. 

One of the other issues of LSE implementation is that the coarse grids should 
capture the flow feature adequately for efficient extrapolation. At the same time these 
coarse grids should not be so fine as to increase the computational cost. As the grid 
resolution is increased the solution reaches an asymptotic trend; i.e., there is no 
noticeable change in the flow feature with further refinement of the grid. Unless the 
coarse grids are properly chosen, the performance of LSE will not be optimal. Judicious 
selection of the base grids is an issue that deserves to be investigated. 

Design Optimization 

In this dissertation the RSAs were first generated and then used for the sensitivity 
and trade-off analyses. Hence although the global sensitivity analyses provided incite into 
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the influence of design variables on the individual objective variability, it was not used 
for dimensionality reduction. Additionally the injector design problem has small number 
of design variables and hence there is not a serious requirement for such a reduction. For 
future design studies the following direction is suggested to tackle the curse of 
dimensionality. Whenever possible the preliminary design study can be conducted using 
a reduced order model, including semi-empirical and experimental data. This reduced 
model can be coupled with DOE, RSM and global sensitivity tools to identify the 
influence of design variables on objective variability. This information can then be used 
to reduce the dimensionality of the design space. Higher fidelity models can then be used 
to analyze the designs in this reduced design space. Additionally, the reduced design 
space will help in reducing the number of search directions during the optimization study. 
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